龙格库塔法的c++编程

gghhjj

1. #include<stdlib.h>
   
2. #include<stdio.h>
   
3. /*n表示几等分，n+1表示他输出的个数*/
   
4. int RungeKutta(double y0,double a,double b,int n,double *x,double *y,int style,double (*function)(double,double))
   
5. {
   
6. double h=(b-a)/n,k1,k2,k3,k4;
   
7. int i;
   
8. // x=(double*)malloc((n+1)*sizeof(double));
   
9. // y=(double*)malloc((n+1)*sizeof(double));
   
10. x[0]=a;
    
11. y[0]=y0;
    
12. switch(style)
    
13. {
    
14. case 2:
    
15. for(i=0;i<n;i++)
    
16. {
    
17. x[i+1]=x[i]+h;
    
18. k1=function(x[i],y[i]);
    
19. k2=function(x[i]+h/2,y[i]+h*k1/2);
    
20. y[i+1]=y[i]+h*k2;
    
21. }
    
22. break;
    
23. case 3:
    
24. for(i=0;i<n;i++)
    
25. {
    
26. x[i+1]=x[i]+h;
    
27. k1=function(x[i],y[i]);
    
28. k2=function(x[i]+h/2,y[i]+h*k1/2);
    
29. k3=function(x[i]+h,y[i]-h*k1+2*h*k2);
    
30. y[i+1]=y[i]+h*(k1+4*k2+k3)/6;
    
31. }
    
32. break;
    
33. case 4:
    
34. for(i=0;i<n;i++)
    
35. {
    
36. x[i+1]=x[i]+h;
    
37. k1=function(x[i],y[i]);
    
38. k2=function(x[i]+h/2,y[i]+h*k1/2);
    
39. k3=function(x[i]+h/2,y[i]+h*k2/2);
    
40. k4=function(x[i]+h,y[i]+h*k3);
    
41. y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;
    
42. }
    
43. break;
    
44. default:
    
45. return 0;
    
46. }
    
47. return 1;
    
48. }
    
49. double function(double x,double y)
    
50. {
    
51. return y-2*x/y;
    
52. }
    
53. //例子求y'=y-2*x/y(0<x<1);y0=1;
    
54. /*
    
55. int main()
    
56. {
    
57. double x[6],y[6];
    
58. printf("用二阶龙格-库塔方法\n");
    
59. RungeKutta(1,0,1,5,x,y,2,function);
    
60. for(int i=0;i<6;i++)
    
61. printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]);
    
62. printf("用三阶龙格-库塔方法\n");
    
63. RungeKutta(1,0,1,5,x,y,3,function);
    
64. for(i=0;i<6;i++)
    
65. printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]);
    
66. printf("用四阶龙格-库塔方法\n");
    
67. RungeKutta(1,0,1,5,x,y,4,function);
    
68. for(i=0;i<6;i++)
    
69. printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]);
    
70. return 1;
    

复制代码

assist

谢谢！

lzh316

还可以，简单点，如果中间加一个迭代

风花雪月

原帖由 lzh316 于 2007-1-22 19:40 发表
还可以，简单点，如果中间加一个迭代

你说的中间加一个迭代是什么意思？

风花雪月

原帖由 lzh316 于 2007-1-22 19:40 发表
还可以，简单点，如果中间加一个迭代

一般这类算法已经很少自己写程序了，大部分在算法库中都能够找到

insects

有没有谁会编个二阶用龙格库塔法解非线性微分方程的程序给我看看，用matlab语言的，方程是


d2θe(t)/dt2+(1/τ1+K*(τ2/τ1)*cosθe(t))* dθe(t)/dt+k/τ1*sinθe(t)-△ω0/τ1=0

风花雪月

原帖由 insects 于 2007-5-16 10:25 发表
有没有谁会编个二阶用龙格库塔法解非线性微分方程的程序给我看看，用matlab语言的，方程是


d2θe(t)/dt2+(1/τ1+K*(τ2/τ1)*cosθe(t))* dθe(t)/dt+k/τ1*sinθe(t)-△ω0/τ1=0

把你的方程转换到状态空间就可以了

注意版规：Matlab问题请到Matlab讨论区发帖询问

sssssxxxxx921

有没有谁会用VC编的啊

风花雪月

原帖由 sssssxxxxx921 于 2007-7-10 21:13 发表
有没有谁会用VC编的啊

什么意思？上面的代码稍加修改就能在VC下运行

佩刀

（用C/C++）用龙格库塔法如何编程微分方程组dx/dt=-10x+10y      dy/dt=30x-y-xz     dz/dt=(-8/3)z+xy求大神帮忙

佩刀

dx/dt=-10x+10y dy/dt=30x-y-xz         dz/dt=(-8/3)z+xy