定步长四阶经典公式 解决数值积分

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
[\frac{dy}{dt} = f(t, y)]
# O4 f2 k0 ^0 h2 j4 \7 X1 p2 y& k( v这个方法的迭代公式如下：
% z5 K, k$ o/ B9 n[k1 = h \cdot f(tn, yn)]
[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
$ ]% \* l5 \7 V' A[k3 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k2}{2})]
9 R7 o* `- \# B; X9 N9 i! \1 Q[k4 = h \cdot f(tn + h, yn + k_3)]
[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m
   # O6 Q4 k. M) i\" d' J2 I
2. 
   ( F3 i% g2 i- O- N
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)/ `- r2 J# c( v, y! C4 {/ D
   
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');
   3 B, l6 ]# r8 I6 V  A0 D
5.    disp('function z=f(x,y)');
   # G% G/ T4 _0 z9 N9 ~6 R
6.    disp('z=y-2*x/y;');9 b  F1 r) e+ N2 n  G: W
   
7.    disp('并将该文件保存在work文件夹下');
   - R+ [  e/ K\" s7 O2 I9 |4 D  j7 v
8. end
   8 q4 ^\" u0 S6 g\" h
9. - {, i+ P+ m& @, ^$ b5 h
   
10. X1=input('请输入求解区间的左端点X1=');
    : t( b6 k3 X- N2 X  i) {3 ?5 u1 P
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');( W8 J( {( e7 N3 Q: L# e9 i
    
12. Xn=input('请输入求解区间的右端点Xn=');, L6 v5 h9 \  y& q7 G
    
13. h=input('请输入求解步长h=');! Z$ v\" p; U3 Q
    
14. / S. H+ ?! A1 y1 T5 Z
    
15. X=X1;7 a# L: ?2 J3 S
    
16. Y=Y1;                                                        %运算初始点4 n\" d8 z/ w5 b+ j( h1 G* V9 z
    
17. n=0;                                                         %节点序号变量置零
    : m\" ~3 A& F8 Z; h! s% D) {
18. 
    + |: k5 [7 I% B3 c& q
19. while X<=Xn-h
    6 L/ I4 r# |: E
20.     K1=f(X,Y);
    / P* R4 D4 A; V. [
21.     K2=f(X+h/2,Y+K1*h/2);
    ; A4 V/ h\" c) s( u$ u
22.     K3=f(X+h/2,Y+K2*h/2);
    \" S3 @1 F, G/ G* L- L2 P
23.     K4=f(X+h,Y+K3*h);
    / \  i( Z, o. X8 X; ?
24.     X=X+h;# `4 b5 I4 l: I) W1 ?
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式, X6 Y\" R& c  u- y. z
    
26.     n=n+1;                                                   %节点序号加1- A/ {, C/ n% \8 k
    
27. ( V' z+ B) h( j: N$ k9 R
    
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);# S  e. R( x* w; \' Z' f% e/ r( Y
    
29.     plot(X,Y,'o')
    # u5 ?\" f9 J\" T! e# z
30.     hold on
    : [$ s8 x5 g+ p
31. end

复制代码

1. function z=f(x,y)
   0 H$ D# \\" z( U7 a4 p$ Q: b3 N
2. z=y-2*x/y;

复制代码

3 l% m  l7 v4 }$ J: G/ U( O0 ]