龙格库塔求解微分方程数值解

kkkttt

工程中很多的地方用到龙格库塔求解微分方程的数值解，
龙格库塔是很重要的一种方法，尤其是四阶的，精确度相当的高。

此代码只是演示求一个微分方程
/*y'=y-2x/y,x∈[0,0.6]
/*y(0)=1
的解，要求解其它的微分方程，可以自己定义借口函数，退换程序里面的函数：

float f(float , float)即可；

1. /*
   
2. /*编程者：刘艮平
   
3. /*完成日期：2004.5.29
   
4. /*e.g:y'=y-2x/y,x∈[0,0.6]
   
5. /*　　y(0)=1
   
6. /*使用经典四阶龙格-库塔算法进行高精度求解　
   
7. /* y(i＋1)＝yi＋( K1+ 2*K2 +2*K3+ K4)/6
   
8. /* K1=h*f(xi,yi)
   
9. /* K2=h*f(xi+h/2,yi+K1/2)
   
10. /* K3=h*f(xi+h/2,yi+K2/2)
    
11. /* K4=h*f(xi+h,yi+K3)
    
12. */
    
13. #include <stdio.h>
    
14. #include <math.h>
    
15. 
    
16. float f(float den,float p0)  //要求解的微分方程的右部的函数 e.g: y-2x/y
    
17. {
    
18.         float rus;
    
19. //   den=w/(W0+sl);
    
20. //   rus=k*A*f/Wp*sqrt(RTp)-(k-1)*.....
    
21.         rus=p0-2*den/p0;
    
22.         return(rus);
    
23. }
    
24. 
    
25. //float fden()
    
26. //{
    
27. 
    
28. //}
    
29. 
    
30. 
    
31. void main()
    
32. {
    
33.         float x0; //范围上限
    
34.         int x1;   //范围下限
    
35.         float h;  //步长
    
36.         int n;    //计算出的点的个数
    
37.         float k1,k2,k3,k4;
    
38.     float y[3]; //用于存放计算出的常微分方程数值解
    
39.         int i=0;
    
40.         int j;
    
41. 
    
42. //以下为函数的接口
    
43.         printf("intput x0:");
    
44.     scanf("%f",&x0);
    
45. 
    
46.         printf("input x1:");
    
47.         scanf("%f",&x1);
    
48. 
    
49.     printf("input y[0]:");
    
50.     scanf("%f",&y[0]);
    
51. 
    
52.         printf("input h:");
    
53.     scanf("%f",&h);
    
54. 
    
55. // 以下为核心程序
    
56.     n=((x1-x0)/h);
    
57.         n=3;
    
58.        
    
59.     for(j=0;j<n;j++)
    
60.         {   
    
61. 
    
62.       k1=h*f(x0,y[i]); //求K1
    
63.       k2=h*f((x0+h/2),(y[i]+k1/2)); //求K2
    
64.       k3=h*f((x0+h/2),(y[i]+k2/2)); //求K3
    
65.       k4=h*f((x0+h),(y[i]+k3)); //求K4
    
66. 
    
67.           y[i+1]=y[i]+((k1+2*k2+2*k3+k4)/6); //求yi+1
    
68.           x0+=float(0.2);
    
69.       printf("y[%f]=%f\n",x0,y[i+1]);
    
70.       ++i;
    
71.          
    
72.         }
    
73. 
    
74. }

复制代码

[ 本帖最后由 风花雪月 于 2006-10-11 15:45 编辑 ]

pheigenbau

我想学偏微分方程的数值解法

lihuaifu

谢谢~等我仔细看看有没有关于非线性微分方程组的问题!

yut_o

楼主有没有关于RK解微分方程组的程序啊

风花雪月

原帖由 yut_o 于 2006-10-11 15:20 发表
楼主有没有关于RK解微分方程组的程序啊

这个就是rk法解微分方程的C语言程序