matlab chính xác yếu tố quyết định vấn đề

Question

tôi có chương trình sau đây

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 

giải pháp phân tích của hệ thống trên, và cùng program written in fortran đưa ra giá trị của omegan bằng 16,3818 và 32,7636 (giá trị fortran; phân tích khác nhau một chút, nhưng họ đang ở đâu đó).

Vì vậy, bây giờ tôi tự hỏi ... tôi đang đi sai với điều này? Tại sao MATLAB không đưa ra kết quả mong đợi?

(điều này có lẽ là một cái gì đó khủng khiếp đơn giản, nhưng nó đem lại cho tôi đau đầu)

Answer 1

câu trả lời mới:

Bạn có thể điều tra vấn đề này sử dụng phương trình mang tính biểu tượng, mà mang lại cho tôi câu trả lời đúng:

>> 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 

Tôi nghĩ bạn không chỉ chọn giá trị trong vòng lặp đủ gần với các giải pháp cho omegan o r ngưỡng của bạn về việc yếu tố quyết định gần bằng không quá nghiêm ngặt. Khi tôi cắm vào các giá trị cho đến a, cùng với omegan = 16.3819 (đó là giá trị gần nhất với một giải pháp vòng lặp của bạn sản xuất), tôi có được điều này:

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

ans = 

    2.765476845475786e-005 

Đó là vẫn lớn hơn trong biên độ tuyệt đối so với 1e-10.

Answer 2

Bạn đang tìm quá ít giá trị xác định vì Matlab đang sử dụng hàm xác định khác (hoặc một số lý do khác giống như độ chính xác của dấu phẩy động liên quan đến hai phương pháp khác nhau). Tôi sẽ cho bạn thấy rằng Matlab về cơ bản là cung cấp cho bạn các giá trị chính xác và một cách tốt hơn để tiếp cận vấn đề này nói chung.

Trước tiên, hãy lấy mã của bạn và thay đổi nó một chút.

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)) 

Tất cả những gì tôi đã thực hiện thực sự được đăng nhập các giá trị của các yếu tố quyết định cho tất cả các giá trị của omegan và vẽ các bản ghi của những giá trị yếu tố quyết định là một hàm của omegan. Dưới đây là cốt truyện:

Bạn nhận thấy ba dips lớn trong đồ thị. Hai trùng với kết quả của bạn là 16.3818 và 32.7636, nhưng cũng có thêm một kết quả mà bạn bị thiếu (có thể do điều kiện của yếu tố quyết định nhỏ hơn 1e-10 quá thấp ngay cả khi mã Fortran của bạn nhận được nó). Do đó, Matlab cũng nói với bạn rằng đó là những giá trị của omegan mà bạn đang tìm kiếm, nhưng vì yếu tố quyết định được xác định theo một cách khác nhau trong Matlab, các giá trị không giống nhau - không đáng ngạc nhiên khi xử lý ma trận kém . Ngoài ra, nó có thể đã làm với Fortran bằng cách sử dụng phao chính xác đơn như người khác nói. Tôi sẽ không xem xét lý do tại sao họ không phải vì tôi không muốn lãng phí thời gian của tôi vào đó. Thay vào đó, hãy nhìn vào những gì bạn đang cố gắng làm và thử một cách tiếp cận khác.

Bạn, như tôi chắc chắn rằng bạn đã biết, đang cố gắng để tìm ra giá trị riêng của ma trận

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

, cài đặt chúng bằng

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

và giải quyết cho omegan. Đây là cách tôi đã đi về nó:

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 

này cung cấp cho bạn tất cả bốn giải pháp - không chỉ hai mà bạn đã tìm thấy với mã Fortran của bạn (mặc dù một trong số họ, không, là ngoài phạm vi thử nghiệm của bạn cho omegan). Nếu bạn muốn giải quyết điều này bằng cách kiểm tra yếu tố quyết định trong Matlab, như bạn đã cố gắng thực hiện, thì bạn sẽ phải chơi với giá trị mà bạn đang kiểm tra giá trị tuyệt đối của yếu tố quyết định nhỏ hơn. Tôi đã nhận nó để làm việc cho một giá trị của 1e-4 (nó đã cho 3 giải pháp: 16.382, 28.374, và 32.764).

Xin lỗi vì giải pháp dài như vậy, nhưng hy vọng nó sẽ giúp ích.

Cập nhật:

Trong khối đầu tiên của tôi về mã trên, tôi đã thay thế

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

với

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

đó là thuật toán mà det được cho là sử dụng theo quy định của Matlab tài liệu. Bây giờ, tôi có thể sử dụng 1E-10 làm điều kiện và nó hoạt động. Vì vậy, có lẽ Matlab không tính toán định thức chính xác như các tài liệu nói? Đây là loại đáng lo ngại.

Answer 3

Tôi đặt câu trả lời này là vì tôi không thể dán điều này vào nhận xét: Đây là cách Matlab tính toán determinant. Tôi giả sử các lỗi làm tròn số đến từ tính là sản phẩm của nhiều yếu tố chéo trong U.

Algorithm

Các yếu tố quyết định được tính toán từ yếu tố tam giác thu được bằng cách khử Gauss

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