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

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.

Answer 1

[Torsten]

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

[Hiroki ENDO]

Thank you for your all advices and hospitality.

I am sorry to late to reply.

This code was worked very well.

Thanks again.

Answer 2

[Walter Roberson]

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.

[Hiroki ENDO]

Thank you for your advice.

I am going to try.

[Bruno Luong]

What theory can ensure the order of eigen value using symbolic is "right"?

[Torsten]

I guess the general solution formula for polynomials of degree <= 4 ensures it in this case.

[Bruno Luong]

Yes but not because using the single-expression formula might not imply the output "right" (change smoothly)

Simple example of root of y^3=x

x=linspace(-1,1,1000)

x = 1×1000

-1.0000 -0.9980 -0.9960 -0.9940 -0.9920 -0.9900 -0.9880 -0.9860 -0.9840 -0.9820 -0.9800 -0.9780 -0.9760 -0.9740 -0.9720 -0.9700 -0.9680 -0.9660 -0.9640 -0.9620 -0.9600 -0.9580 -0.9560 -0.9540 -0.9520 -0.9499 -0.9479 -0.9459 -0.9439 -0.9419

yref=x.^(1/3);

hold on
y1=yref; y1(x<0)=y1(x<0).*exp(2i*pi/3);  plot(y1,'.')
y2=yref; y2(x<0)=y2(x<0).*exp(-2i*pi/3);  plot(y2,'x')
y3=yref;  plot(y3,'o')
legend('y1','y2','y3')

The right one is y1 (blue small dot), the one that has close formula changes.

[Torsten]

So you think that the three formulae for the three different branches don't follow the roots continuously ?

[Bruno Luong]

Not differentiable continuously.

I think using consistent fomal radical formula doesn't ensure there is no branch swapping either.

Answer 3

[Torsten]

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

[Hiroki ENDO]

Thank you for update your answer!

That's very helpfull for me. Thanks.

So, i tried to understand your code.

I have a question about the section of %Compute d(lambda_k)/d(A(i,j)).

I want a time derivative of lambda as d(lambda)/dt, not d(lambda)/dA.

I cannot understand that code.

I ran this in MATLAB. But, I had an error; Unrecognized function or variable 'A'.

So, I changed the code of the function like bellow.

Then I ran the Mtalab. But, the calculation never ends.

[T,Av] = ode15i(@(t,Av, Avd)Adotxx(t, Av, Avd, d, A), tspan, tens2vec(Azero), tens2vec(Adotzero));
function res = Adotxx(t, Av, Avd, d, A)

All functions in a script must be closed with an 'end'.

Error in connector.internal.fevalMatlab

Error in connector.internal.fevalJSON

[Torsten]

I have a question about the section of %Compute d(lambda_k)/d(A(i,j)).

I want a time derivative of lambda as d(lambda)/dt, not d(lambda)/dA.

You must calculate d(lambda)/dt as

d(lambda)/dt = sum_i,j d(lambda)/d(A(i,j)) * d(A(i,j))/dt

as done in the code.

I ran this in MATLAB. But, I had an error; Unrecognized function or variable 'A'.

It seems you must redefine A in function Adotxx as

A = sym('A',[3 3])

(see the change in the code above).

Then I ran the Mtalab. But, the calculation never ends.

Check whether ODE15I starts integrating by outputting the variable t in Adotxx.

[Hiroki ENDO]

Thank you again for your response!

I could understand that. Thanks again.

I will try it soon!

[Hiroki ENDO]

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.

Answer 4

[Torsten]

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

[Hiroki ENDO]

Thank you for your advices and hospitality again!

I am going to try this code after understanding.

Thanks.

[Hiroki ENDO]

I also ran this code. Thank you for providing!!

This calculation takes a lot time. (never ends).

Under your condition, How long does it takes until end?

[Torsten]

The numerical code from above takes about 1 second (under Octave).

[Hiroki ENDO]

Oh. OK.

Thank you for replying.

I will try hard more.