定步长四阶经典公式解决数值积分 - 数学建模社区-数学中国

%四阶经典公式,微分方程为f.m
8 c# u7 V3 r- [, W/ a
. f H5 P M5 o' O3 [
if exist('f.m')==0 %在星号处输入文件名(把星号改为文件名)
r% n# V, r _; \9 e8 [( P# R/ P
disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');
( N( b2 y: g* T) T
disp('function z=f(x,y)');: `& k$ ?" L8 Z2 L) M2 n
disp('z=y-2*x/y;');- @2 `' v" }) _, s( q
disp('并将该文件保存在work文件夹下');) b6 J6 R/ x+ m; e
end
3 d; Y2 v) ~8 W2 v9 Z) T
0 r5 S3 s( m1 n( ?& s$ k
X1=input('请输入求解区间的左端点X1=');
8 n8 F. v' b: _& x
Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');
$ a) V( X+ s/ ]. p ^
Xn=input('请输入求解区间的右端点Xn=');
+ f7 x2 x- o0 ~! p0 Y6 o1 Y5 U/ n& Q
h=input('请输入求解步长h=');
* n- D6 U5 T8 H
' {- N0 I& [" n3 I, g: Y5 p
X=X1; G4 O" [/ W9 v% _9 T/ {
Y=Y1; %运算初始点
# L5 n. e* w) t9 W& N. H
n=0; %节点序号变量置零& G9 t8 l6 }7 m" s3 r Q# e$ j
7 r" ~# Y- n6 J9 _
while X<=Xn-h1 }6 g! c7 J% [9 q7 b$ U9 L
K1=f(X,Y);9 X p" L+ }/ y$ P" ]" c3 S
K2=f(X+h/2,Y+K1*h/2);( v* w& |4 n! O5 M: a- K- v) H
K3=f(X+h/2,Y+K2*h/2);
; o( ^; [2 D' B! u( H
K4=f(X+h,Y+K3*h);$ F9 p$ T9 h1 i" Q& P
X=X+h;
" ~/ ]8 t2 G4 k, P3 V0 c4 x
Y=Y+h*(K1+2*K2+2*K3+K4)/6; %四阶标准的龙格-库塔公式% `% c8 e8 N( ^6 s. U' f, P
n=n+1; %节点序号加1) J$ e5 { ]. S0 E3 @
6 X' ?* C1 A2 x, S$ j# y- D9 R( z. b
fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);9 V. A8 u% X- y# @; x- Z
plot(X,Y,'o')# |/ E. ?5 ^7 i/ J- }
hold on
5 x- E6 \6 I% _' x6 `3 p# n
end

复制代码