目录
一、 引言…………………………………… 3
二、 算法描述……………………………… 4
三、 试验及分析…………………………… 7
四、 自己的特色…………………………… 9
五、 总结…………………………………… 12
六、 程序清单……………………………… 13
一、引言
科学技术中常常要求解常微分方程的定解问题。这类问题最简单的形式是本章将要着重考察的一阶方程的初值问题
我们知道，只要函数适当光滑，譬如关于满足Lipschitz条件
。 (9.1.3)
理论上就可以保证初值问题(9.1.1),(9.1.2)的解存在并且唯一。
虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法。
所谓数值解法，就是寻求在一系列离散节点
上的近似值。相邻两个结点的间距称为步长。今后如果不特别说明，总是假定为定数，这时节点为 。
初值问题(9.1.1),(9.1.2)的数值解法有个基本特点，它们都采用取“步进式”，即求解过程顺着节点的次序一步一步地向前推进。描述这类算法，只要给出用已知信息计算的递推公式。
首先，要对方程(9.1.1)离散化，建立求数值解的递推公式。一类是计算时只用到前一点的值，称为单步法。另一类是用到前面点的值，称为步法。其次，要研究公式的局部截断误差和阶，数值解与精确解的误差估计及收敛性，还有递推公式的计算稳定性等问题。
二、算法描述
1、Euler method
我们知道，在平面上，微分方程(1.1)的解称作它的积分曲线。积分曲线上一点的切线斜率等于函数的值，如果按函数在平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致。
基于上述几何解释，我们从初始点出发，先依方向场在该点的方向推进到上一点，然后再从依方向场的方向推进到上一点，循此前进推出一条折线(图9.1)。
一般地，设已做出该折线的顶点，过依方向场的方向再推进到，显然两个顶点，的坐标有关系
即
(9.2.1)
这就是著名的欧拉(Euler)公式。若初值已知，则依公式(9.2.1)可逐步算出
2、Backward Euler
假设，即顶点落在积分曲线上，那么，按欧拉方法做出的折线便是过点的切线，从图形上看，这样定出的顶点明显地偏离了原来的积分曲线，可见欧拉方法是相当粗糙的。
为了分析计算公式的精度，通常可用泰勒展开将在处展开，则有
在的前提下，，于是可得到欧拉法(9.2.1)的公式误差
， (9.2.3)
称为此方法的局部截断误差。
如果对方程(9.1.1)从到积分,得
． (9.2.4)
右端积分用左矩形公式近似，再以代替，代替也得到(9.2.1)，局部截断误差也是(9.2.3)．
如果在(9.2.4)中右端积分用右矩形公式近似，则得另一个公式
， (9.2.5)
称为后退的欧拉法．
3、Crank-Nicholson
为得到比欧拉法精度高的计算公式，在等式(9.2.4)右端积分中若用梯形求积公式近似，并用代替，代替，则得
， (9.2.7)
称为梯形方法．
梯形方法是隐式单步法，可用迭代法求解．同后退的欧拉方法一样，仍用欧拉方法提供迭代初值，则梯形法的迭代公式为
(9.2.8)
为了分析迭代过程的收敛性，将(9.2.7)式与(9.2.8)相减，得
，
于是有
，
式中为对满足Lipschitz常数，如果选取充分小，使得
，
则当时有，这说明迭代过程(9.2.8)是收敛的．
三、实验及分析
测试例子：设置长h=0.1,步数n=10。
运行环境：VC++6.0
1、Euler method
2、Backward Euler
3、Crank-Nicholson
四、自己的特色
1、对算法经行了多个区间的比较并且分析误差原因、
通过上面三个区间的比较之后我们发现在算法的区间在增大以后，其误差会变的很大，经过我们组内的讨论，发现其函数的原函数为一个指数函数为 exp(1.0*t*t*t/3.0);所以在增大区间的时候会出现很大的误差。
2、对精度进行研究
这是通过三种平台进行的测试从而很好的为我们的编写程序需要的精度经行了估计，其中我们不能发现其中有很多值是有区别的，比如long double的值在不同的平台上是不一样的这是因为。ANSI C标准规定了double变量存储为 IEEE 64 位（8 个字节）浮点数值，但并未规定long double的确切精度。所以对于不同平台可能有不同的实现。这样在VC6.0上面的测试数据为第一个，可以看到最精确的应该为double。
测试代码如下
#include<iostream>
#include<string>
#include <limits>
using namespace std;
int main()
{
cout << "type: \t\t" << "************size**************"<< endl;
cout << "bool: \t\t" << "所占字节数：" << sizeof(bool);
cout << "\t最大值：" << (numeric_limits<bool>::max)();
cout << "\t\t最小值：" << (numeric_limits<bool>::min)() << endl;
………………
return 0;
}
3、改进了Euler公式
显示欧拉公式计算工作量小，但精度低。向后欧拉法、Crank-Nicholson方法虽提高了精度，但为隐式公式，需用迭代法求解，计算工作量大。综合两者便可改进Euler公式。 先用Euler公式求出一个初步的近似值，称为预测值，它的精度不高，再用Crank-Nicholson公式对它校正一次，即迭代一次，求得，称为校正值。即得到改进的Euler公式：
可以证明该公式的精度是二阶，这是一种显示方式
实例：
为改进：
改进后：
从误差的分析看来截距误差减小了
五、总结
这次大家通过交流学到了很多，我将从以下几点经行说明
1、大家对算法经行了细致的研究
经过上三次的编写程序之后我们对各自的编程的特点和基本思路有了一个清晰的了解，大家在交流的时候也是明显的速度加快了许多，同时对微分的认识也是加强了很多。对于这次的编程我们组提出了也尽量做到时间和空间也精确读的最优化。于是大家在课下不仅仅是对课本上的思想在程序经行了实现，而且在实现之后大家相互的交流之后，对算法的时间空间上面进行了优化，在这其中大家试着去做到如何精确，做了许多超出课本上要就的测试，得出为什么越大误差会很大，并且探求出是函数的原因。其次，有的同学提出了对经行两次运算从而更加精确的确定其微分的值。大家除了课本从算法的优化和对积分的认识上面有了很大的提高。
2、大家积极的参与到后期的报告活动中
在这次中，大家表现出前所未有的热情，争先恐后的去做报告PPT等任务，力求做的最好。在即上三次的分工之后，我们这次决定再次相互的交换一下工作。我们最终的程序或许没有解决我们最初讨论出的问题，但是我想大家参与了就是一种最大的收获。这次的分工是：陈俊桦、韩学武、肖大军和梅旭由四位同学负责修改程序，王艳琴和向胜男经行程序的测试和算法的描述工作，彭亚琼、刘宇、方正写word文档，张钟文和杨耀鹏负责制作PPT ，肖大军和梅旭写台上的细化工作。
六、程序清单
#include<stdio.h>
#include<math.h>
double f(double t,double y)
{
return t*t*y;
}
double fr(double t)
{
return exp(1.0*t*t*t/3.0);
}
int EulerMethod()
{
int i;
double Y,t=0.0,y=1.0,dx=0.1,Error[15];//从（0,1）点开始
printf("\n");
printf("***************************Euler Method***************************\n");
printf("No.intervals\tEuler method\t Error\t Error Ratio\n");
for(i=1; i<=10; i++)
{
t=(i-1)*dx;
y=y+dx*f(t,y);//微分近似值
Y=fr(t);//真实值
Error[i]=Y-y;
if(i==1)
{
printf("%d\t %.10f\t %.10f\t \n",i,y,Error[i]);
}
if(i>1)
{
printf("%d\t %.10f\t %.10f\t%.10f \n",i,y,Error[i],Error[i]/Error[i-1]);
}
}
printf("\n\n");
return 0;
}
int BackwardEuler()
{
int i;
double Y,t=0.0,y=1.0,dx=0.1,Error[15];
printf("\n");
printf("***********************Backward Eule**********************\n");
printf("No.intervals\tBackward Eule\t Error\t Error Ratio\n");
for(i=1; i<=10; i++)
{
t=i*dx;
y=y+dx*f(t,y);
Y=fr(t);
Error[i]=Y-y;
if(i==1)
{
printf("%d\t %.10f\t %.10f\t \n",i,y,Error[i]);
}
if(i>1)
{
printf("%d\t %.10f\t %.10f\t%.10f \n",i,y,Error[i],Error[i]/Error[i-1]);
}
}
printf("\n\n");
return 0;
}
int GrankNicholson()
{
int i;
double Y,t=0.0,y=1.0,dx=0.1,Error[15],tem;
printf("\n");
printf("*********************Backward Eule**********************\n");
printf("No.intervals\tGrank-Nicholson\t Error\t Error Ratio\n");
for(i=1; i<=10; i++)
{
t=i*dx;
tem=y;
y=y+dx*f(t,y);
y=tem+(dx*f(t,y)+dx*f((i-1)*dx,tem))/2.0;
Y=fr(t);
Error[i]=Y-y;
if(i==1)
{
printf("%d\t %.10f\t %.10f\t \n",i,y,Error[i]);
}
if(i>1)
{
printf("%d\t %.10f\t %.10f\t%.10f \n",i,y,Error[i],Error[i]/Error[i-1]);
}
}
printf("\n\n");
return 0;
}
int main()
{
while(1)
{
int n;
printf("************************************\n");
printf("\t1.Euler method\n");
printf("\t2.Backward Euler\n");
printf("\t3.Grank-Nicholson\n");
printf("\t0. exit\n");
printf("************************************\n");
printf("Please enter your choice(1,2,3,0):\n");
scanf("%d",&n);
switch(n)
{
case 1:
EulerMethod();
break;
case 2:
BackwardEuler();
break;
case 3:
GrankNicholson();
break;
case 0:
return 0;
}
}
return 0;
}
#include<stdio.h> #include<math.h> #define zero 1e-5 double f(double x) { return exp(-x); } int main() { double a,b; double x1,x2; printf("迭代法的实现\n"); printf("例式：x = e(-x)"); printf("请输入区间【】："); scanf("%lf %lf",&a,&b); x1=b; x2=a; while(fabs(x1-x2)>=zero) { x1 = x2; printf("%.7lf ",x1); x2 = f(x1); printf("%.7lf\n",x2); } return 0; }
#include<iostream> #include<cmath> #include<iomanip> using namespace std; double f(double x) //函数f(x)，容易修改 { double y; //x=0.001*x; y=cos(x)-0.5; return y; } double ft(double x) //函数f的导数 { double y; y=(cos(x+0.001)-cos(x))/0.001; return y; } double root(double x) //求解方程 { double a[10],e[10]; int i=0,j; for(j=0;j<10;j++) { a[j]=0; e[j]=0; } while(fabs(f(x))>0.0001) { a[i++]=x; x=x-f(x)/ft(x); } for(j=0;j<i;j++) e[j]=fabs(x-a[j]); for(j=0;j<i;j++) cout<<e[j]<<" "<<e[j+1]/e[j]<<" "<<e[j+1]/(e[j]*e[j])<<endl; return x; } int main() { cout<<"请输入猜想值："; double x,y; cin>>x; y=root(x); cout<<"方程的根为:"<<y<<endl; }
#include<stdio.h> #include<math.h> #define PI 3.14159265 #define zero 0.0001 double f(double x) { return 2*sin(PI*x) + cos(PI*x); } int main() { double a,b; double x1,x2,x3,y; printf("二分法的实现\n"); printf("例子式：2*sin（π*x）+cos（π*x）=0\n"); printf("请输入根所在区间【】:"); scanf("%lf %lf",&a,&b); printf("ak bk xk\n"); x1=a; x2=b; x3 = (x1+x2)/2; while(fabs(x2-x1)>=zero) { printf("%lf %lf %lf\n",x1,x2,x3); if(f(x1)*f(x3)>0 )x1=x3;//x2=x2|| f(x1)*f(b)<0 else if(f(x1)*f(x3)<0)x2=x3;//|| f(x1)*f(a)<0 else printf("\t\t%.6lf\n",x1); x3 = (x1+x2)/2; } printf("\t\t %.6lf\n",x3); return 0; }