Problem facing in solving xdot=Ax system when A is unstable.
11 views (last 30 days)
Show older comments
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.
3 Comments
Accepted Answer
Sam Chak
on 26 Feb 2024
Edited: Sam Chak
on 26 Feb 2024
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)
%% 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)
function dxdt = odefun(t, x)
A = [0 1;
1 0];
dxdt = A*x;
end
3 Comments
Sam Chak
on 27 Feb 2024
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
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)
%% 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
More Answers (0)
See Also
Categories
Find more on Matrix Indexing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!