How to solve differential equation including derivative of eigenvalue of tensor?

26 次查看（过去 30 天）
Hiroki ENDO
Hiroki ENDO 2022-9-24
评论： Bruno Luong 2022-10-28
I am trying to solve this equation by ODE45.
here, A is sencond-order tensor, and is time derivative of eigenvalue of A.
So, I try to wirte like bellow.
tspan=[0:0.01:1];
Azero= [0.4 0.1 0; 0.1 0.4 0.1; 0 0.1 0.2]
[T,Av] = ode45('Adotxx', tspan, tens2vec(Azero)); % Azero converted to vector style, this is another function
%%
function Aderiv = Adotxx(t, Av)
A = vec2tens(Av); % Azero converted to tensor style from vector style
% making eigenvector
[e, lambda] = eig(A);
% must be sorted eigenvectors and eigenvalues
[lambda_index, SortO] = sort(diag(lambda),'descend');
lambda_temp = zeros(3);
e_temp = zeros(3);
for i =1:3
lambda_temp(i,i) = lambda_index(i);
e_temp(1:3,i) = e(1:3,SortO(i));
end
% time derivative of lambda
lambda_dot = diff(lambda);
X = lambda_dot
% differential equation
Aderiv = A + X
end
However, you know, the lambda_dot is not true time derivative. it is just derivative in matrix at one moment like Differentiation - MATLAB & Simulink - MathWorks. lambda_dot become matrix of (2,3) by diff(lambda).
How can I correct the code?
Thank you for your hospitality.

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2022-9-28
编辑：Torsten 2022-9-28
Ran in:
A third variante of the code also seems to work under MATLAB ( for high values of the integration tolerances :-) )
syms s
A = sym('A',[3 3]);
% Define characteristic polynomial
determinant = det(A-s*eye(3));
% Compute eigenvalues
lambda = solve(determinant==0,s,'MaxDegree',3);
tspan = [0:0.01:1];
Azero = [0.4 0.1 0; 0.1 0.4 0.1; 0 0.1 0.2];
L = double(subs(lambda,A,Azero)).';
Azero = [Azero;L];
Azero = reshape(Azero.',[12,1]);
M = [eye(9),zeros(9,3);zeros(3,12)];
alpha = 0.2;
M(1,10) = -alpha;
M(5,11) = -alpha;
M(9,12) = -alpha;
options = odeset('Mass',M,'RelTol',1e-3,'AbsTol',1e-3);
[T,Av] = ode15s(@(t,A)Adotxx(t,A,lambda), tspan, Azero, options); % Azero converted to vector style, this is another function
plot(T,Av(:,1))
Warning: Imaginary parts of complex X and/or Y arguments ignored.
function Aderiv = Adotxx(t, A, lambda)
A = reshape(A,[3 4]).';
Av = A(1:3,:);
L = A(4,:);
A = sym('A',[3 3]);
LL = double(subs(lambda,A,Av)).';
% differential equation
Aderiv(1:3,1:3) = Av;
Aderiv(4,1:3) = L - LL;
Aderiv = reshape(Aderiv.',[12 1]);
end

更多回答（3 个）

Walter Roberson
Walter Roberson 2022-9-24
You need to construct the formula for the eigenvalues of the derivative based on the equation for A. As you have a 3x3 matrix that will possibly involve the roots of a cubic equation. You must write them out in explicit form. This all must be calculated ahead of time.
Then inside you substitute particular numeric values into the equations inside your ode function. It is important that you do not use eig and sort, because those do not use consistent orders and so you would fail the requirements for continuous derivatives of your equations.
  6 个评论
显示 4更早的评论
Torsten
Torsten 2022-10-28
So you think that the three formulae for the three different branches don't follow the roots continuously ?
Bruno Luong
Bruno Luong 2022-10-28
Not differentiable continuously.
I think using consistent fomal radical formula doesn't ensure there is no branch swapping either.

请先登录，再进行评论。

Torsten
Torsten 2022-9-26
编辑：Torsten 2022-9-27
See if it works. I cannot test it.
syms s
A = sym('A',[3 3]);
% Define characteristic polynomial
determinant = det(A-s*eye(3));
% Compute eigenvalues
lambda = solve(determinant==0,s,'MaxDegree',3)
% Compute d(lambda_k)/d(A(i,j))
for i = 1:3
for j = 1:3
for k = 1:3
d(k,i,j) = diff(lambda(k),A(i,j))
end
end
end
tspan = [0:0.01:1];
Azero = [0.4 0.1 0; 0.1 0.4 0.1; 0 0.1 0.2];
Adotzero = zeros(3,3);
[T,Av] = ode15i(@(t,Av, Avd)Adotxx(t, Av, Avd, d), tspan, tens2vec(Azero), tens2vec(Adotzero)); % Azero converted to vector style, this is another function
function res = Adotxx(t, Av, Avd, d)
A = sym('A',[3 3]);
Av = vec2tens(Av);
Avd = vec2tens(Avd);
d_num = zeros(3,3,3);
for i = 1:3
for j = 1:3
for k = 1:3
d_num(k,i,j) = double(subs(d(k,i,j),A,Av));
end
end
end
lambda_dot = zeros(3,1);
for k =1:3
lambda_dot(k) = 0.0;
for i = 1:3
for j = 1:3
lambda_dot(k) = lambda_dot(k) + d_num(k,i,j)*Avd(i,j);
end
end
end
alpha = 1;
f = @(A) A;
res = Avd - f(Av) - alpha*diag(lambda_dot);
res = tens2vec(res);
end
  4 个评论
显示 2更早的评论
Hiroki ENDO
Hiroki ENDO 2022-9-28
I tried that code chaged by your advices.
And I put bellow code in function for derivative. It can see time record of calculation.
t=t
It works. But values of t become smaller and smaller. It never reaches over 0.01.
So, now I am using the profiler of MATLAB. I don't know if it works. I am seeking any ideas.
If you have any ideas, could you please teach me?
Thank you.

请先登录，再进行评论。

Torsten
Torsten 2022-9-27
编辑：Torsten 2022-9-27
Here is a numerical code to solve your problem.
It's risky to use this approach since - as Walter Roberson noted - the ordering of the eigenvalues coming from "eig" might change during the computation. This will cause wrong results.
Thus the suggested symbolic approach from above is more safe.
tspan=[0:0.01:1];
Azero= [0.4 0.1 0; 0.1 0.4 0.1; 0 0.1 0.2];
[~, lambda] = eig(Azero);
lambda=[lambda(1,1),lambda(2,2),lambda(3,3)];
Azero = [Azero;lambda];
for i=1:4
for j=1:3
Av((i-1)*3+j) = Azero(i,j);
end
end
M = [eye(9),zeros(9,3);zeros(3,12)];
alpha = 0.2;
M(1,10) = -alpha;
M(5,11) = -alpha;
M(9,12) = -alpha;
options = odeset('Mass',M);%,'RelTol',1e-8,'AbsTol',1e-8);
[T,Av] = ode15s(@Adotxx, tspan, Av, options); % Azero converted to vector style, this is another function
plot(T,Av(:,12))
%%
function Ad = Adotxx(t, A)
for i=1:3
for j=1:3
Av(i,j) = A((i-1)*3+j);
end
L(i) = A(9+i);
end
[~, lambda] = eig(Av);
lambda = [lambda(1,1),lambda(2,2),lambda(3,3)];
% differential equation
Aderiv(1:3,1:3) = Av;
Aderiv(4,1:3) = L - lambda;
for i=1:4
for j=1:3
Ad((i-1)*3+j) = Aderiv(i,j);
end
end
Ad=Ad.';
end
  4 个评论
显示 2更早的评论

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Introduction to Installation and Licensing 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by