Problem facing in solving xdot=Ax system when A is unstable.
11 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Hemanta Hazarika
am 26 Feb. 2024
Kommentiert: Sam Chak
am 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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1628198/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1628213/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1628218/image.png)
3 Kommentare
Akzeptierte Antwort
Sam Chak
am 26 Feb. 2024
Bearbeitet: Sam Chak
am 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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1628328/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1628333/image.png)
%% 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 Kommentare
Sam Chak
am 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
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1628863/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1628868/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1628873/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1628878/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1628883/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1628888/image.png)
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
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Matrix Indexing finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!