transponer - problema determinante de la precisión de matlab

symbolic determinant matlab

Question

Tengo el siguiente programa

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

La solución analítica del sistema anterior, y el mismo programa escrito en fortran da valores de omegan igual a 16.3818 y 32.7636 (valores fortánicos, analíticos difieren un poco, pero están ahí en alguna parte).

Entonces, ahora me pregunto ... ¿dónde me estoy equivocando con esto? ¿Por qué Matlab no está dando los resultados esperados?

(Esto es probablemente algo terriblemente simple, pero me está dando dolores de cabeza)

Answer 1

Está buscando valores determinantes demasiado pequeños porque Matlab está usando una función determinante diferente (o alguna otra razón como algo relacionado con la precisión del punto flotante involucrado en los dos métodos diferentes). Te mostraré que Matlab básicamente te está dando los valores correctos y una mejor manera de abordar este problema en general.

Primero, tomemos su código y lo modifiquemos ligeramente.

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

Todo lo que he hecho realmente es registrar los valores del determinante para todos los valores de omegan y trazar el registro de esos valores determinantes como una función de omegan. Aquí está la trama:

Usted nota tres caídas importantes en el gráfico. Dos coinciden con los resultados de 16.3818 y 32.7636, pero también hay uno adicional que te faltaba (probablemente porque tu condición del determinante era menos de 1e-10 era demasiado baja incluso para que tu código Fortran la recogiera). Por lo tanto, Matlab también te está diciendo que esos son los valores de omegan que estabas buscando, pero debido a que el determinante se determinó de manera diferente en Matlab, los valores no fueron los mismos, lo que no es sorprendente cuando se trata de matrices mal condicionadas. . Además, probablemente tenga que ver con Fortran usando flotadores de precisión simple como dijo otra persona. No voy a ver por qué no lo están porque no quiero perder el tiempo con eso. En cambio, veamos lo que estás tratando de hacer y prueba un enfoque diferente.

Usted, como estoy seguro de que sabe, está tratando de encontrar los valores propios de la matriz

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

, configúrelos igual a

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

y resuelve para omegan Así es como lo hice:

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

Esto le da las cuatro soluciones, no solo las dos que había encontrado con su código Fortran (aunque una de ellas, cero, estaba fuera de su rango de prueba para omegan). Si quieres resolver esto al verificar el determinante en Matlab, como has estado tratando de hacer, entonces deberás jugar con el valor de que estás comprobando que el valor absoluto del determinante sea menor que. Lo hice funcionar por un valor de 1e-4 (dio 3 soluciones: 16.382, 28.374 y 32.764).

Perdón por una solución tan larga, pero espero que ayude.

Actualizar:

En mi primer bloque de código anterior, reemplacé

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

con

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

que es el algoritmo que supuestamente se usa de acuerdo con los documentos de Matlab. Ahora, puedo usar 1E-10 como condición y funciona. Entonces, ¿quizás Matlab no está calculando el determinante exactamente como dicen los doctores? Esto es un poco molesto

Answer 2

Puse esto como respuesta porque no puedo pegar esto en un comentario: Así es como Matlab calcula el determinante . Supongo que los errores de redondeo provienen del cálculo del producto de múltiples elementos diagonales en U.

Algoritmo

El determinante se calcula a partir de los factores triangulares obtenidos por eliminación gaussiana

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

Answer 3

Nueva respuesta:

Puedes investigar este problema usando ecuaciones simbólicas, que me dan las respuestas correctas:

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

Creo que o bien no estaba eligiendo los valores en su ciclo lo suficientemente cerca de las soluciones para omegan o su umbral de cuán cerca está el determinante de cero es demasiado estricto. Cuando omegan = 16.3819 los valores dados a a , junto con omegan = 16.3819 (que es el valor más cercano a una solución que produce su bucle), obtengo esto:

>> det(subs(a,{omegan,G,J,Jm,d},{16.3819,8e7,77,10540,140/3}))
ans =
2.765476845475786e-005

Que aún es más grande en amplitud absoluta que 1e-10.