本文介紹兩種入門(mén)級(jí)求解微分方程的方法 —— 梯形法與歐拉法。

將上述方程組改寫(xiě)成matlab語(yǔ)言：

function F = fun(t,Y)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%??????????????????????????????????????????????????????????? %
%???程序作者：Miracle ????????（matlab愛(ài)好者公眾號(hào)）???????????%
%                                                            %
%              歡迎關(guān)注matlab愛(ài)好者公眾號(hào)                      %
%                                                            %
% 任何人都可以免費(fèi)無(wú)條件獲取本程序，切勿將本程序用于商業(yè)用途。%
% 程序版權(quán)歸matlab愛(ài)好者公眾號(hào)所有。%
%                                                            %
% 敬告：切勿刪改本聲明部分，否則將自動(dòng)失去本程序的使用權(quán)               %
%                                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 把初值傳給T、T、V、C
T = Y(1);
Tx = Y(2);
V = Y(3);
C = Y(4);
% 設(shè)置對(duì)應(yīng)的微分方程
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?向后歐拉公式

? ?歐拉法求解源代碼? ?

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                            %
%   程序作者：Miracle         （matlab愛(ài)好者公眾號(hào)）           %
%                                                            %
%              歡迎關(guān)注matlab愛(ài)好者公眾號(hào)                      %
%                                                            %
% 任何人都可以免費(fèi)無(wú)條件獲取本程序，切勿將本程序用于商業(yè)用途。%
% 程序版權(quán)歸matlab愛(ài)好者公眾號(hào)所有。%
%                                                            %
% 敬告：切勿刪改本聲明部分，否則將自動(dòng)失去本程序的使用權(quán)               %
%                                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear;clc;close?all;
Delta?=?0.001;???????%定義步長(zhǎng)
t = 0:Delta:50;      %定義自變量t
n = length(t);       %自變量長(zhǎng)度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');

??歐拉法結(jié)果圖??

二、梯形法

? ?梯形法求解源代碼???

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                            %
%   程序作者：Miracle         （matlab愛(ài)好者公眾號(hào)）           %
%                                                            %
%              歡迎關(guān)注matlab愛(ài)好者公眾號(hào)                      %
%                                                            %
% 任何人都可以免費(fèi)無(wú)條件獲取本程序，切勿將本程序用于商業(yè)用途。%
% 程序版權(quán)歸matlab愛(ài)好者公眾號(hào)所有。%
%                                                            %
% 敬告：切勿刪改本聲明部分，否則將自動(dòng)失去本程序的使用權(quán)               %
%                                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear;clc;close?all;
Delta = 0.001;        % 定義步長(zhǎng)
t = 0:Delta:50;       % 定義自變量t
n = length(t);        % 自變量長(zhǎng)度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');

???梯形法結(jié)果圖? ?

感謝Miracle向公眾號(hào)投稿！歡迎更多愛(ài)好、喜歡matlab編程的朋友來(lái)稿，在公眾號(hào)回復(fù)“投稿”了解投稿詳情。

參考資料：

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

圖片來(lái)源：由?Gerd Altmann 在Pixabay上發(fā)布文章來(lái)源地址http://www.zghlxwxcb.cn/news/detail-553223.html