|
|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
MATLAB源程序代码分享:MATLAB实现改进的欧拉法求解常微分方程组! N3 K. i0 e1 O6 }- x, `
function F=f(t,Y)7 O5 N2 v% F5 F% a
% 定义待求的微分方程组
, L5 }. t% \6 t. E. c. u# ]x=Y(1);
2 }' @+ w) m9 M" D' e0 e* C ly=Y(2);1 @, o- f- h; M5 r: `5 g
f1=3*y;" s/ z' U+ ]7 g; a7 g& J
f2=(1-x^2)*y-x;
* U6 y5 J! N4 j w/ _F=[f1;f2];
% P2 R5 W" x& l, s: X) h5 Z7 _: J& iend$ V: \8 M) a4 {6 W: a- B4 o
% B- X1 D; L; O) s3 o& K
3 z( t2 o8 R5 X, S+ k3 u
%% 定义计算的步长, 设置变量的初始值
" S( w: X4 r- C) a# G, ~% d: Yclear;clc;close all3 A" Z+ S/ y8 d% W& u7 ^" ^3 X4 n
Delta=0.001; % 定义步长1 _# u1 c& g5 F: }
t=0: Delta:20; % 定义自变量 t
- T. g. A/ T8 {5 h6 O- U& \n=length(t);
7 W3 l& V8 ^' @' C# Z$ f+ N/ EY(:,1)=[2;0]; % 定义 x y 的初始值
5 c$ O+ k5 j. M/ ^8 b7 H* Q D* L& s6 g/ B! W; n, I% u6 t
%% 自定义改进的欧拉法, 求解微分方程组
- J/ b |/ K3 ^8 a8 _( G' t. Dfor k=1:n-11 E' p3 Z- j4 `% o
F1=f(t(k),Y( :,k));2 G0 {7 [' T3 _: x5 K* Q( R
F2=f(t(k+1),Y( :,k)+Delta*f(t(k),Y( :,k)));
/ o' |& n1 |; i$ h Y(:,k+1)=Y( :,k)+Delta*(F1+F2)/2;3 Z% n4 H& g- Q* A# _6 M- E
end
& H$ k& C3 L% S; x" N1 Yx=Y(1,: );8 j9 X% ]2 d4 e" h, _! Q( u
y=Y(2,: );3 }( c) C: u6 a
, e7 Y8 h) u/ B$ D& m
%% 绘制 x y 的求解结果; a: o- {0 L. W9 A. e5 M& x
figure- O% K! R5 Q$ {2 V' R% V" w
set(gcf,'units','normalized','position',[0.15 0.2 0.7 0.6]); % 设置 figure 窗口的位置和尺寸
& Y9 S: d& M& L8 ]: Wsubplot(1,2,1)
8 f5 V. k3 V: k/ y9 q+ hplot(t,x,'b')
( P' i4 S+ i! T" z( Zxlabel('t') T A2 m2 H- H; x- e) X/ n2 d6 K& X0 F
ylabel('x')5 o' f' u( v8 ]2 p$ _
5 q+ [. a" S9 V4 a$ T, r$ L- V
subplot(1,2,2)
( w- s/ U4 y( Wplot(t,y,'r')
: x; l/ n0 P) [( y+ Bxlabel('t')' @& m3 v& N5 C) I( V$ |4 l: x, d
ylabel('y')
" t# T) A. T {1 X9 M% X: ?( g
+ {- c9 q/ ~+ m9 |& N
|
|