Problem facing in solving xdot=Ax system when A is unstable.

When solving a problem using ode45 for the system A is special unstable matrix. I encountered a situation where the solution x turned out to be on the order of $10^29$. Since I'm interested in the error between two states, theoretically the error should be zero. However, up to 4 decimal places, the values of the states seems to be the same, i.e., for example, and and but e=x1−x2 is not equal to zero. I understand that x1 and x2 have mismatched in decimal points, and the values are amplified by the exponent. How can I restrict these numerical instabilities? Please help.

[Paul]

You're much more likely to get help if you post or attach your A-matrix (can attach a .mat file using the paperclip icon on the Insert menu) and post your code. If the theoretical error is zero, what is the goal of using ode45?

[Hemanta Hazarika]

function xdot=odefcn(t,x,D,A,B,K,N)
A1=kron(eye(N,N),A);
A2=kron(D,B*K);
xdot=(A1-A2)*x;

end

N=3;

This statement is not inside any function.
(It follows the END that terminates the definition of the function "odefcn".)

A=[0 1 0;0 0 1;3 5 6];
%[v,d]=eig(A)
B=[0;0;1];
D=[1 -1 0;-1 2 -1;0 -1 1];
K=[107,71,20]
x_0=rand(3*N,1)
%x_0=rand(6,1);
tspan=0:.01:10;
%opts = odeset('AbsTol',1e-8,'RelTol',1e-4, 'MaxStep',0.01);
e1=eig(A-1*B*K)
e2=eig(A-3*B*K)
[t,x]=ode45(@(t,x) odefcn(t,x,D,A,B,K,N),tspan,x_0);
e112=(x(:,1)-x(:,4))
plot(t,e112)

% In this example the erroe between x1,x4,x7 should be zero for any intial
% condition theoritically I dont able to fix why matlab thrown error to be
% high 

[Hemanta Hazarika]

In this case expected solution x(:,1) and x(:,4) are need to be exactly same theoritically but i am not able get from that is my concern @Paul

Answer 1

[Sam Chak]

Hi @Hemanta Hazarika

This is an educated guess. If the unstable state matrix is and the initial values are , then the error between the two states is perfectly zero.

%% Analytical approach
syms x(t) y(t)
eqns    = [diff(x,t) == y, 
           diff(y,t) == x];
cond    = [x(0)==1, y(0)==1];
S       = dsolve(eqns, cond)

S = struct with fields:
    y: exp(t)
    x: exp(t)

%% Numerical approach
[t, x]  = ode45(@odefun, [0 10], [1 1]);
plot(t, x), grid on,
xlabel('t'), ylabel('\bf{x}')
legend({'$x(t)$', '$\dot{x}(t)$'}, 'interpreter', 'latex', 'fontsize', 16, 'location', 'NW')

%% error between two states
error   = x(:,1) - x(:,2);
norm(error)

ans = 0

function dxdt = odefun(t, x)
    A    = [0 1; 
            1 0];
    dxdt = A*x;
end

[Hemanta Hazarika]

function xdot=odefcn(t,x,D,A,B,K,N)

A1=kron(eye(N,N),A);

A2=kron(D,B*K);

xdot=(A1-A2)*x;

end

N=3;

A=[0 1 0;0 0 1;3 5 6];

%[v,d]=eig(A)

B=[0;0;1];

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

K=[107,71,20]

x_0=rand(3*N,1)

%x_0=rand(6,1);

tspan=0:.01:10;

%opts = odeset('AbsTol',1e-8,'RelTol',1e-4, 'MaxStep',0.01);

e1=eig(A-1*B*K)

e2=eig(A-3*B*K)

[t,x]=ode45(@(t,x) odefcn(t,x,D,A,B,K,N),tspan,x_0);

e112=(x(:,1)-x(:,4))

plot(t,e112)

% In this example the erroe between x1,x4,x7 should be zero for any intial

% condition theoritically I dont able to fix why matlab thrown error to be

% high

[Hemanta Hazarika]

In this case expected solution x(:,1) and x(:,4) are need to be exactly same theoritically but i am not able get from that is my concern @Sam Chak

[Sam Chak]

Hi @Hemanta Hazarika,

Without the analytical solution, how can you be certain that states and are precisely identical, as depicted in the example I provided in my previous answer? Furthermore, considering the random initial values, if differs from , it is impossible for to be equivalent to .

N   = 3;
A   = [0 1 0;0 0 1;3 5 6];
% [v, d] = eig(A)
B   = [0; 0; 1];
D   = [1 -1 0;-1 2 -1;0 -1 1];
K   = [107, 71, 20];

A1  = kron(eye(N,N),A);
A2  = kron(D,B*K);
Aa  = A1 - A2

Aa = 9×9

0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 -104 -66 -14 107 71 20 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 107 71 20 -211 -137 -34 107 71 20 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 107 71 20 -104 -66 -14

eig(Aa);

x_0     = rand(3*N, 1);
%x_0=rand(6,1);
tspan   = 0:1e-6:1.;
% opts = odeset('AbsTol',1e-8,'RelTol',1e-4, 'MaxStep',0.01);
% e1=eig(A-1*B*K)
% e2=eig(A-3*B*K)
[t, x]  = ode45(@(t,x) odefcn(t,x,D,A,B,K,N), tspan, x_0);
plot(t, x(:,1), t, x(:,4)), grid on
legend('x_1', 'x_4', 'fontsize', 16, 'location', 'NW')

e112    = x(:,1) - x(:,4);
norm(e112)

ans = 168.6748

%% Error between x1 and x4
plot(t, e112), grid on, title('Error')

function xdot=odefcn(t,x,D,A,B,K,N)
    A1  = kron(eye(N,N),A);
    A2  = kron(D,B*K);
    xdot= (A1 - A2)*x;
end