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

%四阶经典公式,微分方程为f.m' R2 H& V6 Y& b2 H! h
/ z! G* l, L) O& F& ~ p
if exist('f.m')==0 %在星号处输入文件名(把星号改为文件名)
8 V+ F! s. ]$ i6 b8 R: `. t# i$ q
disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');) {. [. F) x# L
disp('function z=f(x,y)');
4 C4 a$ o1 ]4 k6 ], p# H# V
disp('z=y-2*x/y;');
* F& \- l' a5 W! X
disp('并将该文件保存在work文件夹下');
% t6 G/ D! y9 y. f1 ~+ g! l* e, _
end ' b" ?+ _# d# Q
: D& Q/ q, B$ S
X1=input('请输入求解区间的左端点X1=');
5 G8 y& |! ?/ t' p$ q$ u5 T
Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');
! O1 H) T B7 \+ w* n+ ^& y' t
Xn=input('请输入求解区间的右端点Xn=');. b. ?" x$ q; A2 R" g F8 X7 R. u7 x
h=input('请输入求解步长h=');9 ?6 H! S" q) ]% h" I1 k# f
. k) O } Y. ~
X=X1;
0 j, S$ C: L; N& o( m
Y=Y1; %运算初始点
! @8 N; u- |" c" T& p
n=0; %节点序号变量置零
1 z$ z3 `4 j$ E& Z
) }9 ~+ n+ K0 q6 q- l
while X<=Xn-h
" [/ _" M& S, b5 f
K1=f(X,Y);
; A) b% a) A% y) y
K2=f(X+h/2,Y+K1*h/2);- o- v' O0 y( ^7 P a0 D5 f
K3=f(X+h/2,Y+K2*h/2);
4 {% K, c( b; P% U
K4=f(X+h,Y+K3*h);
4 H1 {/ {. v+ t5 G5 h+ d
X=X+h;
9 @2 g9 ?: t4 T" N; z* H3 s& j
Y=Y+h*(K1+2*K2+2*K3+K4)/6; %四阶标准的龙格-库塔公式
: z# {9 _0 Z7 S( A7 \
n=n+1; %节点序号加1: R3 o+ \! G( i( \* i0 u
" m$ U' f0 C& ?) l! Y& F
fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);4 ?/ m; g0 q( k# M. ~% T& y# A
plot(X,Y,'o')
S2 X' R' v& Q8 s* w1 x8 q
hold on
6 C5 _0 J( D3 _. X; n, i
end

复制代码