|
|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
本帖最后由 House 于 2020-3-12 18:11 编辑 2 t7 S) x3 `- U9 Z+ B U% Y
1 z' n' |# \' Q, `3 y" ]! b) m! q
MATLAB源程序代码分享:MATLAB实现四阶龙格库塔法求解常微分方程组
7 ?2 |) b: w' l- Ffunction F=f(t,Y), M# ^0 t' b# O7 t# ]
% 定义待求的微分方程组
' y8 @& F9 G9 U! B+ Ax=Y(1);
l# |- U$ M+ y! v$ R4 ]y=Y(2);( ?- b/ f9 B* S5 G
f1=-x^2*(5-y);
+ H5 I1 e8 |$ c) }5 S+ \3 uf2=y*(4-3*x);0 j2 O# Y3 R+ \1 x& ~6 k% G1 \( C0 `
F=[f1;f2];
. ]: b; L1 Q4 Y$ fend
. }' x3 B0 l* O. m4 j* O
, O% \, j$ R9 j* b4 f d: q" Z
+ R& f& N0 z; Z9 z) J! h+ R9 Y5 [7 m4 t5 m9 P
1 u6 x7 j/ x) N2 `& D% C2 f3 O& P%% 定义计算的步长, 设置变量的初始值4 w3 {, w: e, y; z# @1 T
clear;clc;close all
4 u: K2 F* a8 E' x" ?Delta=0.001; % 步长: M% c& U9 L/ i3 p5 l& D
t=0: Delta :8; % 定义自变量 t- `8 U4 p: t( F+ Y# S2 W
n=length(t); & h0 G e, }+ A0 s1 c
Y(:,1)=[0.5;3]; % 变量 x y 的初始值* T9 o7 q9 g0 m' l, [8 M) @
8 | H% Q& V9 U" G' ?% E7 [
%% 自定义龙格库塔法, 求解微分方程组* j! Q5 S1 U" ?4 w, I
for k=1: n-1
1 u5 @6 e- |" ^! C% m3 f0 { z1=f(t(k),Y(:,k));8 x$ `. t7 J+ {4 S7 N+ I8 R
z2=f(t(k)+Delta/2,Y(:,k)+z1*Delta/2);
5 X- P7 Y6 h b/ b0 p) m6 ] z3=f(t(k)+Delta/2,Y(:,k)+z2*Delta/2);: ?- f1 T- l; c H: ~
z4=f(t(k)+Delta,Y(:,k)+z3*Delta);5 C$ H2 W, E: K# m1 A4 e }- m* C
Y(:,k+1)=Y(:,k)+Delta*(z1+2*z2+2*z3+z4)/6;
3 N# w* k7 I2 {& R. P% Oend
6 M) E3 l% {/ N+ m2 H+ U0 o/ Wx=Y(1,: );; t V2 ^0 |7 W: B( ]9 b; Y! Z% o
y=Y(2,: ); D1 y8 a5 S# F- U# W0 s
/ Z8 v* ^5 ?: z% Q5 ?%% 绘制 x y 的求解结果
1 j; L7 n9 \/ {/ H' N8 e( Pfigure
1 H" v' x5 h5 ]set(gcf,'units','normalized','position',[0.15 0.2 0.7 0.6]); % 设置 figure 窗口的位置和尺寸
4 G" d! l4 D! I$ b# _9 @; Msubplot(1,2,1)
' N' A. ?3 ]8 Rplot(t,x,'b')" K$ Y* |6 u" W, H2 H2 H" j
xlabel('t'), j) P% U2 d' W4 h" v! z( ~
ylabel('x')' x! q- ~. P6 v# S
4 A" W6 U! D' y! p- n' ^: G
subplot(1,2,2)
& K4 [! u/ r3 i% S7 @" d* bplot(t,y,'r')9 _$ S7 @) s; n6 E6 r# F
xlabel('t')
, w4 R* f; Q* }9 Xylabel('y'); {" z4 K& l; V6 T9 Q0 T' A k- u3 v
1 I0 \8 n; s& Y K1 G
|
|