*题目描述:*求常微分方程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==0
if (fabs(x - low) < 1e-6) //浮点数之间不能直接判断是否相等,则可以用 fabs(f1-f2)<=给定差值精度 来判断f1和f2是否相等
return 1; //f(low)== f(0)==1
else 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;
else
return 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;
}
输出结果: