目录
一、研究问题
二、C++代码
三、计算结果
一、研究问题
前面几节内容介绍了常微分方程有限差分格式的推导。为加强对本专栏知识的理解,从本节开始,我们补充一些具体算例及相应的编程。
欧拉法的原理及推导请参考:
常微分方程算法之欧拉法(Euler)_欧拉公式求解微分方程-CSDN博客https://blog.csdn.net/L_peanut/article/details/137268437?spm=1001.2014.3001.5501 研究问题为
取步长为0.1。
二、C++代码
#include <cmath>
#include <stdlib.h>
#include <stdio.h>
int main(int argc, char *argv[])
{
int i,N;
double a,b,h,y0,err,maxerr;
double *x,*y;
double f(double x, double y);
double exact(double x);
a=0.0; //求解区域的左端点
b=1.0; //求解区域的右端点
N=10; //总的剖分数
h=(b-a)/N; //步长
//动态分配长度为(N+1)的数组,存放节点坐标
x=(double *)malloc(sizeof(double)*(N+1));
for(i=0;i<=N;i++)
x[i]=a+i*h;
//动态分配长度为(N+1)的数组,存放对应节点的数值解
y=(double*)malloc(sizeof(double)*(N+1));
y0=1.0; //初值
y[0]=y0; //初值
maxerr=0.0;
for(i=0;i<N;i++)
{
y[i+1]=y[i]+h*f(x[i],y[i]); //欧拉法
err=fabs(y[i+1]-exact(x[i+1])); //计算各节点处误差
printf("x[%d]=%.4f, y[%d]=%f, exact=%f, err=%f.\n",i+1,x[i+1],i+1,y[i+1],exact(x[i+1]),err);
if(err>maxerr)
maxerr=err;
}
printf("The max error is %f.\n",maxerr); //打印最大误差
return 0;
}
//右端项函数
double f(double x, double y)
{
return y-2*x/y;
}
//精确解
double exact(double x)
{
return sqrt(1.0+2*x);
}
三、计算结果
x[1]=0.1000, y[1]=1.100000, exact=1.095445, err=0.004555.
x[2]=0.2000, y[2]=1.191818, exact=1.183216, err=0.008602.
x[3]=0.3000, y[3]=1.277438, exact=1.264911, err=0.012527.
x[4]=0.4000, y[4]=1.358213, exact=1.341641, err=0.016572.
x[5]=0.5000, y[5]=1.435133, exact=1.414214, err=0.020919.
x[6]=0.6000, y[6]=1.508966, exact=1.483240, err=0.025727.
x[7]=0.7000, y[7]=1.580338, exact=1.549193, err=0.031145.
x[8]=0.8000, y[8]=1.649783, exact=1.612452, err=0.037332.
x[9]=0.9000, y[9]=1.717779, exact=1.673320, err=0.044459.
x[10]=1.0000, y[10]=1.784771, exact=1.732051, err=0.052720.
The max error is 0.052720.
从计算结果中可看出,欧拉法的计算误差还是比较大的。