前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >欧拉法与梯形法求解微分方程【含matlab源代码】

欧拉法与梯形法求解微分方程【含matlab源代码】

作者头像
巴山学长
发布2021-07-30 10:38:38
3.1K0
发布2021-07-30 10:38:38
举报
文章被收录于专栏:巴山学长
本文介绍两种入门级求解微分方程的方法 —— 梯形法与欧拉法。

将上述方程组改写成matlab语言:

代码语言:javascript
复制
function F = fun(t,Y)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                            %
%   程序作者:Miracle         (matlab爱好者公众号)           %
%                                                            %
%              欢迎关注matlab爱好者公众号                      %
%                                                            %
% 任何人都可以免费无条件获取本程序,切勿将本程序用于商业用途。%
% 程序版权归matlab爱好者公众号所有。%
%                                                            %
% 敬告:切勿删改本声明部分,否则将自动失去本程序的使用权               %
%                                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 把初值传给T、T、V、C
T = Y(1);
Tx = Y(2);
V = Y(3);
C = Y(4);
% 设置对应的微分方程
f1 = 1 - 0.1*T + 0.5*T*(1 - (T + Tx)/1000) - 0.0014*V*T;
f2 = 0.0014*V*T - 0.9*Tx - 0.03*Tx*C;
f3 = 3.09375*Tx - (3+0.007*T)*V;
f4 = 0.03*Tx*C-0.06*C;
% 放在一起
F = [f1;f2;f3;f4];
end

一、欧拉法

1.1 向前欧拉公式

1.2 向后欧拉公式

欧拉法求解源代码

代码语言:javascript
复制

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                            %
%   程序作者:Miracle         (matlab爱好者公众号)           %
%                                                            %
%              欢迎关注matlab爱好者公众号                      %
%                                                            %
% 任何人都可以免费无条件获取本程序,切勿将本程序用于商业用途。%
% 程序版权归matlab爱好者公众号所有。%
%                                                            %
% 敬告:切勿删改本声明部分,否则将自动失去本程序的使用权               %
%                                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;clc;close all;
Delta = 0.001;       %定义步长
t = 0:Delta:50;      %定义自变量t
n = length(t);       %自变量长度n
Y(:,1) = [20.7172;2;3.1478*10^5;122.1667]; %定义T T* V C的初始值

%% 自定义欧拉法,求解微分方程组
for k = 1:n-1
    %向前欧拉法
    %Y(:,k+1) = Y(:,k) + Delta*f(t(k),Y(:,k));
    Y(:,k+1) = Y(:,k) + Delta*f(t(k),Y(:,k));
    Y(:,k+1) = Y(:,k) + Delta*f(t(k+1),Y(:,k+1));
end
% 给T、T*、V、C赋值
T = Y(1,:);
T_xing = Y(2,:);
V =  Y(3,:);
C =  Y(4,:);
%% 绘制图像
figure;
set(gcf,'units','normalized','position',[0.15 0.2 0.7 0.6]);
subplot(2,2,1);
plot(t,T,'linewidth',1);xlabel('t');ylabel('T');title('T');
subplot(2,2,2);
plot(t,T_xing,'linewidth',1);xlabel('t');ylabel('T^*');title('T^*');
subplot(2,2,3);
plot(t,V,'linewidth',1);xlabel('t');ylabel('V');title('V');
subplot(2,2,4)
plot(t,C,'linewidth',1);xlabel('t');ylabel('C');title('C');

欧拉法结果图

二、梯形法

梯形法求解源代码

代码语言:javascript
复制

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                            %
%   程序作者:Miracle         (matlab爱好者公众号)           %
%                                                            %
%              欢迎关注matlab爱好者公众号                      %
%                                                            %
% 任何人都可以免费无条件获取本程序,切勿将本程序用于商业用途。%
% 程序版权归matlab爱好者公众号所有。%
%                                                            %
% 敬告:切勿删改本声明部分,否则将自动失去本程序的使用权               %
%                                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;clc;close all;
Delta = 0.001;        % 定义步长
t = 0:Delta:50;       % 定义自变量t
n = length(t);        % 自变量长度n
Y(:,1) = [20.7172;2;3.1478*10^5;122.1667];%定义T T* V C的初始值

%% 自定义梯形公式法,求解微分方程组
for k = 1:n-1    
    Y(:,k+1) = Y(:,k) + Delta*f(t(k),Y(:,k));
    Y(:,k+1) = Y(:,k) + Delta*(f(t(k),Y(:,k))+f(t(k+1),Y(:,k+1)));
end
% 给T、T*、V、C赋值
T = Y(1,:);
T_xing = Y(2,:);
V =  Y(3,:);
C =  Y(4,:);
%% 绘制图像
figure;
set(gcf,'units','normalized','position',[0.15 0.2 0.7 0.6]);
subplot(2,2,1);plot(t,T,'linewidth',1);xlabel('t');ylabel('T');title('T');
subplot(2,2,2);plot(t,T_xing,'linewidth',1);xlabel('t');ylabel('T^*');title('T^*');
subplot(2,2,3);plot(t,V,'linewidth',1);xlabel('t');ylabel('V');title('V');
subplot(2,2,4);plot(t,C,'linewidth',1);xlabel('t');ylabel('C');title('C');

梯形法结果图

感谢Miracle向公众号投稿!欢迎更多爱好、喜欢matlab编程的朋友来稿,在公众号回复“投稿”了解投稿详情。

参考资料:

[1] https://blog.csdn.net/weixin_42141390/article/details/110184743 [2] https://blog.csdn.net/misskissC/article/details/8913941

图片来源:由 Gerd Altmann 在Pixabay上发布

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2021-07-27,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 巴山学长 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档