matlab練習程式（線性常微分方程組矩陣解）

Dsp Tian發表於2024-05-16
Matlab矩陣

之前有透過ode和simulink解線性常微分方程組。

除了上面兩種方法，線性常微分方程組還可以透過矩陣的方法求解。

比如下面這個之前使用的方程組：

x'' = x' - x + y' -z'

y'' = y' - y - x'

z'' = z' - z + x'

可以寫成下面矩陣形式：

設這個矩陣為A，那麼解可以表示為如下形式：

可以直接透過matlab的expm函式求解。

或者對A做特徵值分解[e,v]=eig(A)，然後得到特徵值的exp對角陣，最後透過尤拉公式求解。

matlab程式碼如下：

clear all;close all;clc;
X0=[3;1;-1;1;0.5;0];

fvdp = @(T,Y) [Y(4);        %x'
    Y(5);                   %y'
    Y(6);                   %z'
    Y(4)-Y(1)+Y(5)-Y(6);    %x''
    Y(5)-Y(2)-Y(4);         %y''
    Y(6)-Y(3)+Y(4)];        %z'']

[T,Y] = ode45(fvdp, [0:0.1:5],X0);

for i=1:6
    plot(T,Y(:,i),'r-o');
    hold on;
end

A = [0 0 0 1 0 0;
    0 0 0 0 1 0;
    0 0 0 0 0 1;
    -1 0 0 1 1 -1;
    0 -1 0 -1 1 0;
    0 0 -1 1 0 1];

[e,v] = eig(A);

re=[];
for t = 0:0.1:5

    M = e*(exp(real(v)*t).*(cos(imag(v)*t)+1i*sin(imag(v)*t)))*inv(e);
    %M = e*expm(v*t)*inv(e);
    %M = expm(A*t);
    Y = M*X0;

    re=[re;Y'];
end

Y = re-re(1,:);

for i=1:6
    plot(T,real(Y(:,i))+X0(i),'b*');
end
grid on;

結果和ode解法是一致的：

相關文章