c语言-四阶龙格-库塔法
1 #include<stdio.h> 2 #include<math.h> 3 #define n 14 4 //double func1(double x, double y); 5 double func2(double x, double y); 6 int main(){ 7 double a = 0.0, b = 1.4;//求解区间为[a,b] 8 double h, m = 0.0;//h为步长，取0.1； m为y的初值。 9 double x[n + 1] = { 0.0 }, y[n + 1] = { 0.0 };//初始化数组 10 h = (b - a) / n;//求步长 11 x[0] = a; 12 y[0] = m; 13 printf("要求解的常微分方程为：dy/dx=1+y^2,初值为y(0)=0.0;\n"); 14 printf("下面为预估校正法结果：\n"); 15 //下面为预估校正法（即向前欧拉法和梯形法的结合，也称改进的欧拉法） 16 for (int i = 1; i <= n; i++) 17 { 18 y[i] = y[i - 1] + h*(1 + y[i - 1] * y[i - 1]); 19 x[i] = a + i*h; 20 y[i] = y[i - 1] + h*0.5*((1 + y[i - 1] * y[i - 1]) + (1 + y[i] * y[i])); 21 printf("x[%d]=%e , y[%d]=%e\n", i, x[i], i, y[i]); 22 } 23 printf("下面为四阶龙格-库塔法结果：\n"); 24 double k1 = 0.0, k2 = 0.0, k3 = 0.0, k4 = 0.0; 25 for (int i = 1; i <= n; i++) 26 { 27 x[i] = 0.0, y[i] = 0.0; 28 } 29 x[0] = a; 30 y[0] = m; 31 for (int i = 1; i <= n; i++) 32 { 33 k1 = h*func2(x[i - 1], y[i - 1]); 34 k2 = h*func2(x[i - 1] + h / 2.0, y[i - 1] + 0.5*k1); 35 k3 = h*func2(x[i - 1] + h / 2.0, y[i - 1] + 0.5*k2); 36 k4 = h*func2(x[i - 1] + h, y[i - 1] + k3); 37 y[i] = y[i - 1] + (1.0 / 6.0)*(k1 + 2 * k2 + 2 * k3 + k4); 38 printf("x[%d]=%e , y[%d]=%e\n", i, x[i], i, y[i]); 39 x[i] = a + i*h; 40 } 41 return 0; 42 } 43 /* 44 double func1(double x, double y)//dy/dx=x-y+1,func1=x-y+1; 45 { 46 double f = 0; 47 f = x - y + 1; 48 return f; 49 } 50 */ 51 double func2(double x, double y)//dy/dx=1+y*y,func2=1+y*y; 52 { 53 double f = 0.0; 54 f = 1 + y*y; 55 return f; 56 }
-*--------------------------
R-K法一般形式：
