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

上述系统的解析解,并且在同一写入Fortran程序给出到16.3818和32.7636 omegan相等的值(FORTRAN值;分析有点不同,但他们在那里的地方).

那么,现在我想知道......我在哪里错了？为什么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。