Problem facing in solving xdot=Ax system when A is unstable.
11 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Hemanta Hazarika
il 26 Feb 2024
Commentato: Sam Chak
il 27 Feb 2024
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 Commenti
Risposta accettata
Sam Chak
il 26 Feb 2024
Modificato: Sam Chak
il 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 Commenti
Sam Chak
il 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
Più risposte (0)
Vedere anche
Categorie
Scopri di più su Matrix Indexing in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!