*题目描述：*求常微分方程y’=y-2x/y && y(0)=1的数值解。解析解为：y=sqrt(1+2x)
编译环境vc++6.0：
代码如下：
/*1、y'=y-2x/y && y(0)=1;2、解析解为y=sqrt(1+2x)*/#include <stdio.h>#include <math.h>#define N 5//离散点个数,分别取5,8,10，对应步长0.2,0.125,0.1#define low 0//定义域x区间下限#define high 1.0//定义域x区间上限#define h (high-low)/N //步长/*************************4、解析解实现**********************/double AnalyticalSolution(double x) {return sqrt(1 + 2 * x);}/*************************1、显示欧拉实现**********************/double EulerMethod(double x){//若x==0if (fabs(x - low) < 1e-6) //浮点数之间不能直接判断是否相等，则可以用 fabs(f1-f2)<=给定差值精度 来判断f1和f2是否相等return 1;//f(low)== f(0)==1else return EulerMethod(x - h) + h * (EulerMethod(x - h) - (2 * (x - h) / EulerMethod(x - h)));}/*************************2、隐式欧拉实现**********************///隐式欧拉实现由下面两个函数组成double ImplicitEulerMethod_1(double x)/*隐示欧拉第一次迭代,用显示欧拉计算f(x)的初值，带入隐式欧拉方程 f(x)=f(x-1)+h*(f(x)-2*x/f(x))*/{if (fabs(x - low) < 1e-6)//浮点数之间不能直接判断是否相等，则可以用 fabs(f1-f2)<=给定差值精度 来判断f1和f2是否相等return 1;elsereturn ImplicitEulerMethod_1(x - h) + h * (EulerMethod(x) - 2 * x / EulerMethod(x)); //迭代}double ImplicitEulerMethod_2(double x)/*隐示欧拉第2--n次迭代，用第一次迭代的结果f(x)作为初值带入隐式欧拉方程进行第二次迭代，第二次迭代的结果f(x)作为初值带入隐式欧拉方程进行第三次迭代.....直到本次迭代的结果y1与上次的迭代结果y0，满足|y1-y0|<1e-6,则认为收敛迭代结束*/{if (fabs(x - low) < 1e-6) //浮点数之间不能直接判断是否相等，则可以用 fabs(f1-f2)<=给定差值精度 来判断f1和f2是否相等return 1;else {double y0 = ImplicitEulerMethod_1(x), y1 = ImplicitEulerMethod_1(x);//将第一次迭代的结果作为初值带入第二次迭代do {y0 = y1;y1 = ImplicitEulerMethod_2(x - h) + h * (y0 - 2 * x / y0);} while (fabs(y1 - y0) < 1e-6);//本次迭代结果减去上一次迭代结果fabs(y1-y0)<=1e-6 认为是收敛return y1;}}/*************************3、两步欧拉实现**********************/double DoubleStepEuler(double x)//两步欧拉{if (fabs(x - low) < 1e-6)//若x==0 第一个初值给定return 1;else if (fabs(x - low - h) < 1e-6) //若x==low+h 第二个初值用显示欧拉求出return EulerMethod(x);else return DoubleStepEuler(x - 2 * h) + 2 * h * (DoubleStepEuler(x - h) - 2 * (x - h) / DoubleStepEuler(x - h));}/*************************四阶龙格-库塔法**********************/double C_R_K_M(double x){if (fabs(x - low) < 1e-6)//若x==0 第一个初值给定return 1;else {double k1 = C_R_K_M(x - h) - 2 * (x - h) / C_R_K_M(x - h);double k2 = C_R_K_M(x - h) + h / 2.0 * k1 - (2 * (x - h) + h) / (C_R_K_M(x - h) + h / 2.0 * k1);double k3 = C_R_K_M(x - h) + h / 2.0 * k2 - (2 * (x - h) + h) / (C_R_K_M(x - h) + h / 2.0 * k2);double k4 = C_R_K_M(x - h) + h * k3 - (2 * (x - h) + h) / (C_R_K_M(x - h) + h * k3);return C_R_K_M(x - h) + h / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4);}}/*************************4、格式化输出函数**********************/void prit(double (*p)(double x)){int count=0;double i;for (i = low; high - i >0 ; i += h) {printf("y(%0.3lf)=%0.5lf\t", i, p(i));count++;if (count % 4 == 0)printf("\n");}printf("y(%0.3lf)=%0.5lf\t", i, p(i));//high==i时printf("\n\n\n");}int main(){printf("解析解格式化输出:\n");prit(AnalyticalSolution);printf("显示欧拉法格式化输出:\n");prit(EulerMethod);printf("隐示欧拉法格式化输出:\n");prit(ImplicitEulerMethod_2);printf("两步欧拉法格式化输出:\n");prit(DoubleStepEuler);printf("四阶龙格-库塔法格式化输出:\n");prit(C_R_K_M);return 0;}
输出结果：
