|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
MATLAB源程序代码分享:MATLAB实现改进的欧拉法求解常微分方程组9 t& ]+ D* W! W
function F=f(t,Y)4 e, A, R* P0 G4 ~! y" _& F
% 定义待求的微分方程组/ q1 h" \* E# n, q+ |. f
x=Y(1);& f3 [$ u+ j& A2 T
y=Y(2);
$ ]1 c+ [: t$ {# w0 Gf1=3*y;
% O2 b% s! ?. s [! b/ hf2=(1-x^2)*y-x;
' S2 @& W( W8 ]! _# V% xF=[f1;f2];) B) Z" W; G0 s+ r
end4 V c4 {/ _1 \- ?9 A* k
2 F' }: v0 f3 v
- }6 W, A$ C4 T9 K
%% 定义计算的步长, 设置变量的初始值( {+ c4 B/ G. j% w1 _
clear;clc;close all
# ^8 r# C& q L$ g8 l/ r S' |Delta=0.001; % 定义步长
2 K0 x9 W! I5 b8 r- D( P. v) z1 `t=0: Delta:20; % 定义自变量 t
' Y" O1 y' o" u6 g. a( M* f( Bn=length(t);
: t$ R! d2 j' ], v& Z4 B4 S5 tY(:,1)=[2;0]; % 定义 x y 的初始值
1 ?' Q8 R! t1 D+ X5 o- i" A( `' W m1 L( X
%% 自定义改进的欧拉法, 求解微分方程组
! Z g5 c" j) R' o, v1 Ofor k=1:n-1
# e1 x- q: p7 h, J/ O8 d4 h F1=f(t(k),Y( :,k));9 h! o4 y0 Q! Z
F2=f(t(k+1),Y( :,k)+Delta*f(t(k),Y( :,k)));9 K4 R* R. v! w5 D
Y(:,k+1)=Y( :,k)+Delta*(F1+F2)/2;
7 w, p6 v3 V9 K6 V- d) F' z) Uend
+ i( @! p6 q7 M2 l3 ?& g; ^- I2 R3 D- Ax=Y(1,: ); I: G( y/ `( R3 s
y=Y(2,: );" X8 K( G6 U0 C4 S7 ]$ J, @
' V+ K" h* k2 q- ~9 k
%% 绘制 x y 的求解结果- z o6 f' D5 t, O* O9 w
figure& k9 z" \! `- o- S6 h* A6 z
set(gcf,'units','normalized','position',[0.15 0.2 0.7 0.6]); % 设置 figure 窗口的位置和尺寸
$ a; |( |! o: y$ @+ V1 R' D; i, Fsubplot(1,2,1)
i V. |! C: A# @( fplot(t,x,'b')
' A9 ]4 J' n9 e* q/ h9 |0 l+ \xlabel('t')3 ~* P/ l' J, e
ylabel('x')) T/ a0 v8 `0 d& {1 c; k
- J+ p/ w& j; _* ]3 ?1 ~
subplot(1,2,2)
+ S) f, i4 T( ^% z0 g9 ]/ b3 D$ \, rplot(t,y,'r')
- R* p( _$ x" o8 gxlabel('t')7 z) S& Z, M8 j" P1 ?9 Z
ylabel('y')
: u0 G/ T) ?3 K4 T, Y" ^, s! b$ r& [9 i" H5 F
|
|