微分方程组求解（Euler法、RK法）

时间：2020-4-4

微分方程组求解

ode45

Euler

4阶R-K

matlab代码实例

% li
% 2020-3-31
clear; clc; close all
%% euler
f = @(t,Y)[(-8/3)*Y(1)+Y(2)*Y(3),(-10)*Y(2)+10*Y(3),-Y(1)*Y(2)+28*Y(2)-Y(3)];
h = 0.001;t = 0:h:100;N = length(t);
X = zeros(3,N);DX = zeros(3,1);
x0 = [0 0 10^(-10)]';X(:,1) = x0;
for i = 1:N-1
    DX = f(t(i),X(:,i));
    X(:,i+1) = X(:,i)+h*DX';
end
figure(1)
plot3(X(1,:),X(2,:),X(3,:))
%% order4 RK
f = @(t,Y)[(-8/3)*Y(1)+Y(2)*Y(3),(-10)*Y(2)+10*Y(3),-Y(1)*Y(2)+28*Y(2)-Y(3)];
h = 0.001;t = 0:h:100;N = length(t);
X = zeros(3,N);
x0 = [0 0 10^(-10)]';X(:,1) = x0;
for i = 1:N-1
    K1 = f(t(i),X(:,i));
    K2 = f(t(i)+h/2,X(:,i)+(h/2)*K1');
    K3 = f(t(i)+h/2,X(:,i)+(h/2)*K2');
    K4 = f(t(i)+h,X(:,i)+h*K3');
    X(:,i+1) = X(:,i) + (h/6)*(K1'+2*K2'+2*K3'+K4');
end
figure(2)
plot3(X(1,:),X(2,:),X(3,:))
%% ode45
[t,x] = ode45(@lorenz,[0 100],x0);
plot3(x(:,1),x(:,2),x(:,3))

对结果的比较

h = 0.01
ode45与euler的差距：

ode45与4阶R-K的差距：
效果上，似乎Euler方法更好一些