求变步长四阶龙格库塔法的c程序！！

ypp207

一直困扰小弟，这里就拜各位大虾了！！:handshake

konggu168

请参阅《C常用算法程序集》第二版的P177；或第三版的P294.如果找不到这两本书，发短消息，或邮件konggu168@sohu.com

风花雪月

1. /* vsrk4.c: Fourth-order Runge-kutta method with variable step. */
   
2. #define N 100
   
3. #include "vsrk4.c"
   
4. double f( double x, double y );
   
5. extern int vsrk4( double hm[2], double x01[2],
   
6.                  double x[N], double y[N], double eps);
   
7. void main( )
   
8. {
   
9. int i,k;
   
10. static double y0,eps,x01[2],hm[2],x[N],y[N];
    
11. eps=1.0e-8;
    
12. x01[0]=0.0;     /* x01[0]=a */
    
13. x01[1]=1.0;     /* x01[1]=b */
    
14. hm[0]=0.0005;   hm[1]=0.1;
    
15. y0=1.0;         y[0]=y0;
    
16. x[0]=x01[0];
    
17. system("cls");
    
18. printf("Variable step Fourth-order Runge-Kutta method:\n");
    
19. k=vsrk4( hm, x01, x, y, eps );
    
20. printf("     xn           yn\n");
    
21. for(i=0;i<=k;i++)
    
22.    printf("%12.8f %12.8f\n",x[i],y[i]);
    
23. getch();
    
24. }
    
25. /********* function f( x, y ) ********/
    
26. double f( double x, double y )
    
27. {
    
28. double f;
    
29. f=y-2.0*x/y;
    
30. return(f);
    
31. }

复制代码

1. /* vsrk4.c: Fourth-order Runge-kutta method with variable step. */
   
2. #include <stdio.h>
   
3. #include <stdlib.h>
   
4. #include <math.h>
   
5. extern double f( double x, double y );
   
6. double rk4( double x, double y, double hh );
   
7. int vsrk4(double hm[2],double x01[2],double x[],double y[],double eps)
   
8. {
   
9. int i;
   
10. double xx,h,hmax,hmin,x0,x1,yn,ynh,delta;
    
11. hmin=hm[0]; hmax=hm[1];
    
12. x0=x01[0];  x1=x01[1];
    
13. h=hmax;
    
14. xx=x0;
    
15. i=0;
    
16. x[0]=xx;
    
17. while( xx < x1 )
    
18.   {
    
19. again:  if( xx+h > x1 ) h=x1-xx;
    
20.   yn=rk4(xx,y[i],h);
    
21.   ynh=rk4(xx,y[i],h/2.0);
    
22.   ynh=rk4(xx+h/2.0,ynh,h/2.0);
    
23.   delta=fabs((yn-ynh)/15.0);
    
24.   if( ( delta < eps )||( h == hmin ) )
    
25.     {
    
26.     xx=xx+h;
    
27.     i=i+1;
    
28.     x[i]=xx;
    
29.     y[i]=ynh;
    
30.     if( i >= 99 ) goto endd;
    
31.     if( delta < 0.05*eps )
    
32.       {
    
33.       h=2.0*h;
    
34.       if( h > hmax ) h=hmax;
    
35.       }
    
36.     }
    
37.   else
    
38.     {
    
39.     h=h/2.0;
    
40.     if( h < hmin ) h=hmin;
    
41.     goto again;
    
42.     }
    
43.   }
    
44. endd:return(i);
    
45. }
    
46. /**************** fourth order Runge-Kutta *************************/
    
47. double rk4( double xn, double yn, double h )
    
48. {
    
49. double ynp,k1,k2,k3,k4;
    
50. k1=h*f(xn,yn);
    
51. k2=h*f(xn+0.5*h,yn+0.5*k1);
    
52. k3=h*f(xn+0.5*h,yn+0.5*k2);
    
53. k4=h*f(xn+h,yn+k3);
    
54. ynp=yn+(k1+2.0*k2+2.0*k3+k4)/6.0;
    
55. return(ynp);
    
56. }

复制代码

煜宸0922

回复 3 # 风花雪月 的帖子

大侠，有变步长龙格库塔法的matlab程序吗？急求！不胜感激

liliangbiao

回复 4 # 煜宸0922 的帖子

Matlab还用编写RK程序吗？软件自带的ode45比起任何已经共享的C程序来说，精确多了，那可是一帮专门搞算法的大侠编写的供调用的function啊

lnicholas

我也来蹭个程序学学。。。

zzf8811

回复 3 # 风花雪月 的帖子

这个代码有误。。。