微分方程组求解(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方法更好一些

©️2020 CSDN 皮肤主题: 我行我“速” 设计师: Amelia_0503 返回首页
实付0元
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值