Problem facing in solving xdot=Ax system when A is unstable.
11 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
Sam Chak
2024-2-26
编辑:Sam Chak
2024-2-26
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 个评论
Sam Chak
2024-2-27
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
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Matrix Indexing 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!