欧拉方法：
#include <stdio.h>
#include <math.h>
#define N 10//离散点个数取10，对应步长0.1
#define low 0//定义域x区间下限
#define high 1.0//定义域x区间上限
#define h (high-low)/N//步长
double AnalyticalSolution(double x)//准确解
{
return pow((1+pow(x,2)),1.0/3);
}
double EulerMethod(double x)//欧拉公式
{
//若x==0
if(fabs(x-low)<1e-6)//浮点数之间不能直接判断是否相等，则可以用 fabs(f1-f2)<=给定差值精度 来判断f1和f2是否相等
return 1;//y(low)==y(0)==1
else
{
return EulerMethod(x-h)+h*(2.0/3*(x- h)*pow(EulerMethod(x-h),-2));
}
}
void prit(double (*p)(double x))//格式化输出
{
int count=0;
double i;
for(i=low+h;high-i>h;i+=h)
{
printf("y(%0.3lf)=%0.6lf ", i, p(i));
count++;
if(count%5==0)
printf("\n");
}
printf("y(%0.3lf)=%0.6lf\t",i,p(i));//high==i时
printf("\n\n");
}
int main()
{
printf("准确解格式化输出:\n");
prit(AnalyticalSolution);
printf("欧拉方法格式化输出:\n");
prit(EulerMethod);
return 0;
}
改进欧拉方法：
#include <stdio.h>
#include <math.h>
#define N 10//离散点个数取10，对应步长0.1
#define low 0.0//定义域x区间下限
#define high 1.0//定义域x区间上限
#define h (high-low)/N//步长
double AnalyticalSolution(double x)//准确解
{
return pow((1+pow(x,2)),1.0/3);
}
double ImplicitEulerMethod(double x)//改进欧拉方法
{
if(fabs(x-low)<1e-6)//浮点数之间不能直接判断是否相等，则可以用 fabs(f1-f2)<=给定差值精度 来判断f1和f2是否相等
return 1;
else
return ImplicitEulerMethod(x-h)+h/2.0*((2.0/3*(x-h)*pow(ImplicitEulerMe thod(x-h),-2))+(2.0/3*x*pow((ImplicitEulerMethod(x-h)+h*(2.0/3*(x-h)*pow(Implic itEulerMethod(x-h),-2))),-2)));
}
void prit(double (*p)(double x))//格式化输出
{
int count=0;
double i;
for(i=low+h;high-i>h;i+=h)
{
printf("y(%0.3lf)=%0.6lf ", i, p(i));
count++;
if(count%5==0)
printf("\n");
}
printf("y(%0.3lf)=%0.6lf\t",i,p(i));//high==i时
printf("\n\n");
}
int main()
{
printf("准确解格式化输出:\n");
prit(AnalyticalSolution);
printf("改进欧拉方法格式化输出:\n");
prit(ImplicitEulerMethod);
return 0;
}
经典R-K法：
#include <stdio.h>
#include <math.h>
#define N 10//离散点个数取10，对应步长0.1
#define low 0//定义域x区间下限
#define high 1.0//定义域x区间上限
#define h (high-low)/N//步长
double AnalyticalSolution(double x)//准确解
{
return pow((1+pow(x,2)),1.0/3);
}
double R_K(double x)//经典R-K法
{
if(fabs(x-low)<1e-6)//若x==0,第一个初值给定
return 1;
else {
double k1=2.0/3*(x-h)*pow(R_K(x-h),-2);
double k2=2.0/3*((x-h)+h/2.0)*pow((R_K(x-h)+h/2.0*k1),-2);
double k3=2.0/3*((x-h)+h/2.0)*pow((R_K(x-h)+h/2.0*k2),-2);
double k4=2.0/3*((x-h)+h)*pow((R_K(x-h)+h*k3),-2);
return R_K(x-h)+h/6.0*(k1+2*k2+2*k3+k4);//返回y_i+1
}
}
void prit(double (*p)(double x))//格式化输出
{
int count=0;
double i;
for(i=low+h;high-i>h;i+=h)
{
printf("y(%0.3lf)=%0.6lf ", i, p(i));
count++;
if(count%5==0)
printf("\n");
}
printf("y(%0.3lf)=%0.6lf\t",i,p(i));//high==i时
printf("\n\n");
}
int main()
{
printf("准确解格式化输出:\n");
prit(AnalyticalSolution);
printf("经典R-K法格式化输出:\n");
prit(R_K);
return 0;
}
