欧拉法、预估校正法(改进的欧拉法)与四阶龙格库塔法求解常微分方程的数值解C++程序
2016-10-20 13:27
363 查看
以y'=x+y,0<x<1,y(0)=1为例,取步长h=0.1,已知精确值为y=-x-1+2e^x,用来进行精度比较。
#include<stdio.h>
using namespace std;
double cor[10000];
double f(double x,double y)//改写函数
{
return x+y;
}
double correctf(double x)//精确解函数
{
return -x-1+2*exp(x);
}
void Euler(double h,double l,double r,double *a,double *b,double tol)//欧拉法
{
double sum=0;
for(int i=1; i<=tol; i++)
{
b[i]=b[i-1]+h*f(a[i-1],b[i-1]);
sum+=fabs(b[i]-cor[i])/cor[i];
}
for(int i=1; i<=tol; i++)
printf("当x=%lf时,近似解为:%lf,准确解为:%lf\n",a[i],b[i],cor[i]);
printf("精度为:%lf\n\n",sum/tol);
}
void improvedEuler(double h,double l,double r,double *a,double *b,double tol)//改进的欧拉法
{
double b1,sum=0;
for(int i=1; i<=tol; i++)
{
b1=b[i-1]+h*f(a[i-1],b[i-1]);
b[i]=b[i-1]+h/2*(f(a[i-1],b[i-1])+f(a[i],b1));
}
for(int i=1; i<=tol; i++)
printf("当x=%lf时,近似解为:%lf,准确解为:%lf\n",a[i],b[i],cor[i]);
printf("精度为:%lf\n\n",sum/tol);
}
void RungeKutta(double h,double l,double r,double *a,double *b,double tol)//四阶龙格库塔法
{
double k1,k2,k3,k4,sum=0;
for(int i=1; i<=tol; i++)
{
k1=f(a[i-1],b[i-1]);
k2=f(a[i-1]+h/2,b[i-1]+h/2*k1);
k3=f(a[i-1]+h/2,b[i-1]+h/2*k2);
k4=f(a[i-1]+h,b[i-1]+h*k3);
b[i]=b[i-1]+h/6*(k1+2*k2+2*k3+k4);
}
for(int i=1; i<=tol; i++)
printf("当x=%lf时,近似解为:%lf,准确解为:%lf\n",a[i],b[i],cor[i]);
printf("精度为:%lf\n\n",sum/tol);
}
int main()
{
double h,a[10000],b[10000],l,r;
memset(a,0,sizeof(a));
memset(b,0,sizeof(b));
memset(cor,0,sizeof(cor));
printf("请输入步长:");
scanf("%lf",&h);
printf("请输入区间下限:");
scanf("%lf",&l);
printf("请输入区间上限:");
scanf("%lf",&r);
printf("请赋予初始值:");
scanf("%lf",&b[0]);
double tol=(r-l)/h;
for(int i=0; i<=tol; i++)
a[i]=l+i*h;
for(int i=1; i<=tol; i++)
cor[i]=correctf(a[i]);
printf("以下为欧拉法求解结果:\n");
Euler(h,l,r,a,b,tol);
printf("以下为改进的欧拉法求解结果:\n");
improvedEuler(h,l,r,a,b,tol);
printf("以下为四阶龙格库塔法求解结果:\n");
RungeKutta(h,l,r,a,b,tol);
return 0;
}
运行得到:
#include<stdio.h>
using namespace std;
double cor[10000];
double f(double x,double y)//改写函数
{
return x+y;
}
double correctf(double x)//精确解函数
{
return -x-1+2*exp(x);
}
void Euler(double h,double l,double r,double *a,double *b,double tol)//欧拉法
{
double sum=0;
for(int i=1; i<=tol; i++)
{
b[i]=b[i-1]+h*f(a[i-1],b[i-1]);
sum+=fabs(b[i]-cor[i])/cor[i];
}
for(int i=1; i<=tol; i++)
printf("当x=%lf时,近似解为:%lf,准确解为:%lf\n",a[i],b[i],cor[i]);
printf("精度为:%lf\n\n",sum/tol);
}
void improvedEuler(double h,double l,double r,double *a,double *b,double tol)//改进的欧拉法
{
double b1,sum=0;
for(int i=1; i<=tol; i++)
{
b1=b[i-1]+h*f(a[i-1],b[i-1]);
b[i]=b[i-1]+h/2*(f(a[i-1],b[i-1])+f(a[i],b1));
}
for(int i=1; i<=tol; i++)
printf("当x=%lf时,近似解为:%lf,准确解为:%lf\n",a[i],b[i],cor[i]);
printf("精度为:%lf\n\n",sum/tol);
}
void RungeKutta(double h,double l,double r,double *a,double *b,double tol)//四阶龙格库塔法
{
double k1,k2,k3,k4,sum=0;
for(int i=1; i<=tol; i++)
{
k1=f(a[i-1],b[i-1]);
k2=f(a[i-1]+h/2,b[i-1]+h/2*k1);
k3=f(a[i-1]+h/2,b[i-1]+h/2*k2);
k4=f(a[i-1]+h,b[i-1]+h*k3);
b[i]=b[i-1]+h/6*(k1+2*k2+2*k3+k4);
}
for(int i=1; i<=tol; i++)
printf("当x=%lf时,近似解为:%lf,准确解为:%lf\n",a[i],b[i],cor[i]);
printf("精度为:%lf\n\n",sum/tol);
}
int main()
{
double h,a[10000],b[10000],l,r;
memset(a,0,sizeof(a));
memset(b,0,sizeof(b));
memset(cor,0,sizeof(cor));
printf("请输入步长:");
scanf("%lf",&h);
printf("请输入区间下限:");
scanf("%lf",&l);
printf("请输入区间上限:");
scanf("%lf",&r);
printf("请赋予初始值:");
scanf("%lf",&b[0]);
double tol=(r-l)/h;
for(int i=0; i<=tol; i++)
a[i]=l+i*h;
for(int i=1; i<=tol; i++)
cor[i]=correctf(a[i]);
printf("以下为欧拉法求解结果:\n");
Euler(h,l,r,a,b,tol);
printf("以下为改进的欧拉法求解结果:\n");
improvedEuler(h,l,r,a,b,tol);
printf("以下为四阶龙格库塔法求解结果:\n");
RungeKutta(h,l,r,a,b,tol);
return 0;
}
运行得到:
相关文章推荐
- 四阶龙格-库塔法求解常微分方程的初值问题-matlab通用程序
- 全区间积分的阿当姆斯预报校正法(常微分方程组的求解)
- 计算方法之改进的欧拉法计算常微分方程
- 微分方程的数值解法——常微分方程——欧拉法与改进欧拉法(2)
- C++实现经典四阶龙格库塔法解一阶微分方程
- 求解矩阵的最简行矩阵的C++程序代码
- 按照高等代数的传统解法编写c++程序实现N(N<256)元的线性方程组的求解
- C++中读取一个txt数值的程序
- 用户输入M,N值,从1至N开始顺序循环数数,每数到M输出该数值,直至全部输出。写出C/C++程序(链表)
- [算法] [常微分方程] [欧拉法 改进欧拉法 经典R-K算法]
- 数值计算一阶常微分方程求解实现
- 四阶龙格-库塔法求解常微分方程的初值问题
- 一个关于c++数值转换的小程序的思考
- 计算方法之改进的欧拉法计算常微分方程
- 求解常微分方程初值问题之改进Euler法:预报-校正公式
- 【c++程序】求解一元二次方程
- 数值分析常用的几个小程序C++实现
- 数值分析 Gauss-Seidel迭代法求解线性方程组 MATLAB程序实现
- 数值作业:改进欧拉法求常微分方程C语言代码
- 梯形法、龙格库塔法求解微分方程 数值解 笔记