matlab परिशुद्धता निर्धारक समस्या

Question

मैं निम्नलिखित कार्यक्रम

format compact; format short g; clear; clc; 
L = 140; J = 77; Jm = 10540; G = 0.8*10^8; d = L/3; 
for i=1:500000 
omegan=1.+0.0001*i; 

a(1,1) = ((omegan^2)*(Jm/(G*J))*d^2)-2; a(1,2) = 2; a(1,3) = 0; a(1,4) = 0; 
a(2,1) = 1; a(2,2) = ((omegan^2)*(Jm/(G*J))*d^2)-2; a(2,3) = 1; a(2,4) = 0; 
a(3,1) = 0; a(3,2) = 1; a(3,3) = ((omegan^2)*(Jm/(G*J))*d^2)-2; a(3,4) = 1; 
a(4,1) = 0; a(4,2) = 0; a(4,3) = 2; a(4,4) = ((omegan^2)*(Jm/(G*J))*d^2)-2; 

if(abs(det(a))<1E-10) sprintf('omegan= %8.3f det= %8.3f',omegan,det(a)) 
end   
end 

ऊपर प्रणाली का विश्लेषणात्मक समाधान, और एक ही program written in fortran 16.3818 और 32.7636 (fortran मूल्यों के omegan बराबर के मूल्यों बाहर देता है, विश्लेषणात्मक, एक छोटे से अलग है, लेकिन वे कहीं कहीं हैं)।

तो, अब मैं सोच रहा हूं ... मैं इसके साथ गलत कहां जा रहा हूं? Matlab अपेक्षित परिणाम क्यों नहीं दे रहा है?

(यह शायद बहुत सरल कुछ है, लेकिन यह मेरे सिर दर्द दे रही है)

Answer 1

न्यू जवाब:

आप प्रतीकात्मक समीकरण का प्रयोग कर इस समस्या की जांच कर सकते हैं, जो मुझे सही जवाब देता है:

>> clear all    %# Clear all existing variables 
>> format long   %# Display more digits of precision 
>> syms Jm d omegan G J %# Your symbolic variables 
>> a = ((Jm*(d*omegan)^2)/(G*J)-2).*eye(4)+... %# Create the matrix a 
     diag([2 1 1],1)+... 
     diag([1 1 2],-1); 
>> solns = solve(det(a),'omegan') %# Solve for where the determinant is 0 

solns = 

           0 
           0 
      (G*J*Jm)^(1/2)/(Jm*d) 
      -(G*J*Jm)^(1/2)/(Jm*d) 
     -(2*(G*J*Jm)^(1/2))/(Jm*d) 
     (2*(G*J*Jm)^(1/2))/(Jm*d) 
    (3^(1/2)*(G*J*Jm)^(1/2))/(Jm*d) 
-(3^(1/2)*(G*J*Jm)^(1/2))/(Jm*d) 

>> solns = subs(solns,{G,J,Jm,d},{8e7,77,10540,140/3}) %# Substitute values 

solns = 

        0 
        0 
    16.381862247021893 
-16.381862247021893 
-32.763724494043785 
    32.763724494043785 
    28.374217734436371 
-28.374217734436371 

मुझे लगता है कि आप या तो omegan ओ के समाधान के लिए पर्याप्त रूप से अपने लूप में मूल्यों का चयन नहीं कर रहे थे निर्धारक को शून्य के करीब कितना करीब है, यह बहुत सख्त है। जब मैं a को दिया मूल्यों में प्लग, omegan = 16.3819 के साथ (जो एक ही समाधान अपने पाश का उत्पादन के लिए निकटतम मान है), मैं इस मिल:

>> det(subs(a,{omegan,G,J,Jm,d},{16.3819,8e7,77,10540,140/3})) 

ans = 

    2.765476845475786e-005 

कौन सा अभी भी 1e-10 से पूर्ण आयाम में बड़ा है।

Answer 2

आप बहुत कम निर्धारक मानों की तलाश में हैं क्योंकि मैटलैब एक अलग निर्धारक फ़ंक्शन का उपयोग कर रहा है (या दो अन्य तरीकों से जुड़े फ्लोटिंग पॉइंट सटीकता के साथ कुछ अन्य कारण)। मैं आपको दिखाऊंगा कि मैटलैब अनिवार्य रूप से आपको सही मूल्य और सामान्य रूप से इस समस्या से संपर्क करने का एक बेहतर तरीका दे रहा है।

सबसे पहले, चलिए अपना कोड लें और इसे थोड़ा बदल दें।

format compact; format short g; clear; clc; 
L = 140; J = 77; Jm = 10540; G = 0.8*10^8; d = L/3; 
vals = zeros(1,500000); 
for i=1:500000 
    omegan=1.+0.0001*i; 

    a(1,1) = ((omegan^2)*(Jm/(G*J))*d^2)-2; a(1,2) = 2; a(1,3) = 0; a(1,4) = 0; 
    a(2,1) = 1; a(2,2) = ((omegan^2)*(Jm/(G*J))*d^2)-2; a(2,3) = 1; a(2,4) = 0; 
    a(3,1) = 0; a(3,2) = 1; a(3,3) = ((omegan^2)*(Jm/(G*J))*d^2)-2; a(3,4) = 1; 
    a(4,1) = 0; a(4,2) = 0; a(4,3) = 2; a(4,4) = ((omegan^2)*(Jm/(G*J))*d^2)-2; 
    vals(i) = abs(det(a)); 
    if(vals(i)<1E-10) 
     sprintf('omegan= %8.3f det= %8.3f',omegan,det(a)) 
    end 
end 
plot(1.+0.0001*(1:500000),log(vals)) 

सभी कि मैं वास्तव में किया है omegan के सभी मानों के लिए निर्धारक के मूल्यों लॉग इन और omegan के एक समारोह के रूप में उन निर्धारक मूल्यों का लॉग साजिश रची है।

आप ग्राफ में तीन प्रमुख डुबकी नोटिस: यहाँ साजिश है। दो 16.3818 और 32.7636 के परिणामों के साथ मेल खाते हैं, लेकिन एक अतिरिक्त भी है जिसे आप याद कर रहे थे (संभवतः क्योंकि आपके फोरट्रान कोड को लेने के लिए निर्धारक की आपकी स्थिति 1e-10 से भी कम थी)। इसलिए, मैटलैब आपको यह भी बता रहा है कि वे ओमेगन के मूल्य हैं जिन्हें आप ढूंढ रहे थे, लेकिन निर्धारक की वजह से मैटलैब में एक अलग तरीके से निर्धारित किया गया था, मूल्य समान नहीं थे - बुरी तरह से वातानुकूलित मैट्रिक्स से निपटने पर आश्चर्य की बात नहीं । इसके अलावा, इसे शायद फोरट्रान के साथ एकल परिशुद्धता फ्लोट का उपयोग करना होगा जैसा कि किसी और ने कहा था। मैं नहीं देख रहा हूं कि वे क्यों नहीं हैं क्योंकि मैं उस पर अपना समय बर्बाद नहीं करना चाहता हूं। इसके बजाए, देखते हैं कि आप क्या करने की कोशिश कर रहे हैं और एक अलग दृष्टिकोण का प्रयास करें।

आप, के रूप में मुझे यकीन है कि आप जानते हैं, मैट्रिक्स

a = [[-2 2 0 0]; [1 -2 1 0]; [0 1 -2 1]; [0 0 2 -2]]; 

eigenvalues ​​के मिल जाए, सेट उन्हें बराबर

-omegan^2*(Jm/(G*J)*d^2) 

और omegan के लिए हल करने के लिए कोशिश कर रहे हैं रहा हूँ। इस तरह मैं इसके बारे में चला गया है:

format compact; format short g; clear; clc; 
L = 140; J = 77; Jm = 10540; G = 0.8*10^8; d = L/3; 
C1 = (Jm/(G*J)*d^2); 
a = [[-2 2 0 0]; [1 -2 1 0]; [0 1 -2 1]; [0,0,2,-2]]; 
myeigs = eig(a); 
myeigs(abs(myeigs) < eps) = 0.0; 
for i=1:4 
    sprintf('omegan= %8.3f', sqrt(-myeigs(i)/C1)) 
end 

यह आप सभी चार समाधान देता है - न सिर्फ दो है कि आप अपने फोरट्रान कोड के साथ मिल गया था (हालांकि उनमें से एक, शून्य, omegan के लिए अपने परीक्षण सीमा के बाहर था)। यदि आप Matlab में निर्धारक की जांच करके इसे हल करने के बारे में जाना चाहते हैं, जैसा कि आप करने का प्रयास कर रहे हैं, तो आपको उस मान के साथ खेलना होगा जिसे आप निर्धारक के पूर्ण मूल्य से कम होने की जांच कर रहे हैं। मुझे इसे 1e-4 के मान के लिए काम करने के लिए मिला (यह 3 समाधान दिए: 16.382, 28.374, और 32.764)।

इतने लंबे समाधान के लिए खेद है, लेकिन उम्मीद है कि यह मदद करता है।

अद्यतन:

ऊपर कोड की मेरी पहली ब्लॉक में, मैं

vals(i) = abs(det(a)); 

प्रतिस्थापित

[L,U] = lu(a); 
s = det(L); 
vals(i) = abs(s*prod(diag(U))); 

साथ जो एल्गोरिथ्म कि det माना जाता है कि Matlab के अनुसार उपयोग कर रहा है डॉक्स। अब, मैं शर्त के रूप में 1E-10 का उपयोग करने में सक्षम हूं और यह काम करता है। तो शायद मैटलैब निर्धारक की गणना नहीं कर रहा है जैसे दस्तावेज़ कहते हैं? यह परेशान है।

Answer 3

मैंने इसे एक उत्तर के रूप में रखा क्योंकि मैं इसे एक टिप्पणी में पेस्ट नहीं कर सकता: यहां बताया गया है कि मैटलैब determinant की गणना कैसे करता है। मुझे लगता है गोलाई त्रुटियों यू

में कई विकर्ण तत्वों के उत्पाद की गणना एल्गोरिथ्म

निर्धारक त्रिकोणीय गाऊसी उन्मूलन

[L,U] = lu(A) s = det(L)   
%# This is always +1 or -1 
det(A) = s*prod(diag(U))