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

%四阶经典公式,微分方程为f.m
, I( A( U5 M2 P9 E
! R1 j1 |/ b8 |6 }6 ?. G
if exist('f.m')==0 %在星号处输入文件名(把星号改为文件名): {& z1 h P* t6 T' f% i' I" t5 ?
disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');$ c: U0 [/ g/ X
disp('function z=f(x,y)');
( N3 ]9 x9 O$ ]" t
disp('z=y-2*x/y;');# I, W) G8 y0 m) {& ]8 p
disp('并将该文件保存在work文件夹下');
/ t& ~: @ m: P; X
end
$ }& o2 Z5 _! K W' a4 X. R
& m3 E% P$ ?* Q5 t$ v
X1=input('请输入求解区间的左端点X1=');5 Z& [: ~; e. Z& [+ t
Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）'); v" x9 _. t0 G$ z
Xn=input('请输入求解区间的右端点Xn=');
1 q+ y# k9 \/ B5 B& V5 i
h=input('请输入求解步长h=');( J. d) S' n8 a; w2 I; p
# ~" D* D+ E M, _( Q/ r
X=X1;
. u: m/ P) x- n3 U( ~
Y=Y1; %运算初始点" Z4 w3 v" X' H1 ^# T7 o E" u$ E
n=0; %节点序号变量置零7 T+ n. _# u: Q0 @. V- C+ z
' j( R1 e+ Y8 `5 V/ l8 c% }
while X<=Xn-h
3 r; V, d& Q+ f( A
K1=f(X,Y);9 M' h/ \! ^4 k6 ^& \
K2=f(X+h/2,Y+K1*h/2);
- L! N' Z) V- ]9 `" p$ u
K3=f(X+h/2,Y+K2*h/2);
- u S- A; I4 K2 [2 R: Z
K4=f(X+h,Y+K3*h);
* z' l! ~4 p; z% x3 o+ g2 c% e
X=X+h;5 i; m- c2 ?5 s
Y=Y+h*(K1+2*K2+2*K3+K4)/6; %四阶标准的龙格-库塔公式- ?0 ?! c& F! k7 m0 J/ ?9 G( e
n=n+1; %节点序号加1
0 d! Y4 U4 ]# ~7 [ x( g6 V4 I9 z
0 W, m* x5 v$ r# U. F5 G; @5 A
fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);
: H4 x1 K+ d) ]7 i/ P
plot(X,Y,'o')
5 Z: e4 X: ]& [. S. x) f
hold on. q2 {( [9 [1 a
end

复制代码