|
|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
MATLAB源程序代码分享:MATLAB实现改进的欧拉法求解常微分方程组
: Y1 v! B% ~& _8 |& F" V6 Bfunction F=f(t,Y)
+ C h5 }) h+ ]# T% 定义待求的微分方程组
) L$ r e% o* @ {9 K, Z( q- l ?x=Y(1);* F" [" U' S. M
y=Y(2);* C# R. L/ g) p- V5 v4 z
f1=3*y;
9 p. s) f7 p4 {' \$ j: hf2=(1-x^2)*y-x;" U! L# e" t: j6 z* I& ]5 H
F=[f1;f2];" i, \$ o! b; p. E, N7 I
end
8 l: m" O0 y0 u7 w o
; u3 d. w( f! ~ p( m0 C- l0 F. s! N7 P# @( F# l4 R/ n
%% 定义计算的步长, 设置变量的初始值! o3 C% j- w' p" [6 T' |3 l1 A
clear;clc;close all3 z8 m; z# R" k( E4 ?" e F w1 @
Delta=0.001; % 定义步长
) o+ X- `1 r; g, Dt=0: Delta:20; % 定义自变量 t# Z# m I$ |$ |& F( I
n=length(t);
2 s' S9 z+ r f& ^* \2 AY(:,1)=[2;0]; % 定义 x y 的初始值
9 V( r: P( ~8 ?
6 w( ~+ m) O5 u [2 }%% 自定义改进的欧拉法, 求解微分方程组" ]$ R9 T: W' f
for k=1:n-1
# W; G" Z+ c# N% ?' i F1=f(t(k),Y( :,k));0 O, X) v* ]9 M. O( ^4 P8 Q
F2=f(t(k+1),Y( :,k)+Delta*f(t(k),Y( :,k)));
; `$ p+ A- q) f# r6 F Y(:,k+1)=Y( :,k)+Delta*(F1+F2)/2;( }6 h2 `6 e2 u. H% ^" D' M0 h5 R
end+ H9 I8 `7 g+ v. j, K
x=Y(1,: );. g7 E z/ y9 r( F( V
y=Y(2,: );
' t% h6 R( E8 [! g
. G8 d2 n9 y( ?3 W%% 绘制 x y 的求解结果+ `$ P6 |6 G d; b. H! B
figure; V% z# Z# u8 P+ c
set(gcf,'units','normalized','position',[0.15 0.2 0.7 0.6]); % 设置 figure 窗口的位置和尺寸
+ y3 \( P$ A0 M+ z: Osubplot(1,2,1)! ?5 F7 z9 l$ S8 a
plot(t,x,'b'). \- }0 A5 W( Z e6 j
xlabel('t')' m$ t4 {: Q% j4 Q2 u2 w2 G
ylabel('x')
" ?1 U) W9 c" K$ U0 W% ^+ _- a: B' u% F- @2 h4 C; I: [2 {
subplot(1,2,2)
; I- O3 i2 w" ^6 o$ Aplot(t,y,'r')# {& j' P- E* f/ V1 U4 [
xlabel('t')( u& q3 c" w* n0 M5 B
ylabel('y')4 V6 J% v3 C/ _# n* e
' P' _. j( Q" @$ K; s% _
|
|