课题10常微分方程初值问题的数值方法.docx
《课题10常微分方程初值问题的数值方法.docx》由会员分享，可在线阅读，更多相关《课题10常微分方程初值问题的数值方法.docx（12页珍藏版）》请在冰点文库上搜索。
课题10常微分方程初值问题的数值方法
课题10.常微分方程初值问题的数值方法
一．问题提出
（1）利用欧拉方法和改进的欧拉方法求解初值问题：
dy/dx=4*x/y-x*y;y（0）=3;其中0二.问题算法、c语言编程和上机运算结果
1.欧拉法
算法：
输入微分方程f（x，y）；
输入积分部数n；
输入初值x0，y0；
输入步长h；
利用k1=f（xn,yn）
y（n+1）=yn+k1*h；
n=0,1,2……
采取不断循环计算；
输出x1，x2，……xn.
程序：
/*微分方程*/
floatf（x,y）
floatx,y;
{
floatz;
z=4*x/y-x*y;
return（z）;
}
/*欧拉法*/
floatEULAR（f）
float（*f）（）;
{
inti,n;
floatx0,y0,x,y,k1,h;
printf（"\n请输入积分步数n:
"）;
scanf（"%d",&n）;
printf（"\n请输入初值x（0）y（0）:
"）;
scanf（"%f%f",&x0,&y0）;
printf（"\n请输入步长h:
"）;
scanf（"%f",&h）;
printf（"\nxy"）;
printf（"\n%f%f",x0,y0）;
for（i=1;i<=n;i++）
{
x=x0+h;
k1=（*f）（x0,y0）;
y=y0+h*k1;
printf（"\n%f%f",x,y）;
x0=x;
y0=y;
}
}
main（）
{
EULAR（f）;
}
结果：
（1）h=0.1时
（2）h=0.2时
（3）h=0.4时
2.改进的欧拉法
算法：
输入微分方程f（x，y）；
输入积分部数n；
输入初值x0，y0；
输入步长h；
利用k1=f（xn,yn）
k2=f（xn+h,yn+h*k1）
y（n+1）=yn+（k1+k2）/2；
n=0,1,2……
采取不断循环计算；
输出x1，x2，……xn.
程序：
/*微分方程*/
floatf（x,y）
floatx,y;
{
floatz;
z=y-2.0*x/y;
return（z）;
}
/*改进欧拉法*/
floatEULAR（f）
float（*f）（）;
{
inti,n;
floatx0,y0,x,y,k1,k2,h;
printf（"\n请输入积分步数n:
"）;
scanf（"%d",&n）;
printf（"\n请输入初值x（0）y（0）:
"）;
scanf（"%f%f",&x0,&y0）;
printf（"\n请输入步长h:
"）;
scanf（"%f",&h）;
printf（"\nxy"）;
printf（"\n%f%f",x0,y0）;
for（i=1;i<=n;i++）
{
x=x0+h;
k1=（*f）（x0,y0）;
k2=（*f）（x,y0+h*k1）;
y=y0+h*（k1+k2）/2.0;
printf（"\n%f%f",x,y）;
x0=x;
y0=y;
}
}
main（）
{
EULAR（f）;
}
结果：
（1）当h=0.1时
（2）当h=0.2时
（3）当h=0.4时
三.结果分析讨论
1.对比欧拉法，改进的欧拉法和精确解的结果可知，改进的欧拉法所得到结果的精度比欧拉法的大，这是因为改进的欧拉法融入了属于隐式公式的梯形公式，它的计算数值解的精度要比欧拉公式好。
2.从结果分析可以发现，无论是欧拉法还是改进的欧拉法，当步长h取得越小时，所得到的解的精度越大，因此减小步长h可以提高精度。
3.利用欧拉法和改进的欧拉法求解常微分方程初值问题是可行的。
附：
（国庆放假时编的程序）课题4迭代格式的比较
王佳琪学号092071
一．各迭代函数编程及初步分析
1.x=（3*x+1）/（x*x）c语言编程为：
#include
main（）
{
floatx0=1.0,x1;inti;
x1=（3*x0+1）/（x0*x0）;
printf（"%1.5f",x1）;
for（i=1;i<=50&&fabs（x1-x0）>1.0E-5;i++）
{
x0=x1;
x1=（3*x0+1）/（x0*x0）;
printf（"%1.5f",x1）;
}
printf（"\nyidaitimes:
%d",i）;
}
即当x0=1.0时其发散，结果为：
说明此种迭代格式不能得到应求的根。
2.x=（x*x*x-1）/3c语言编程为：
#include
main（）
{
floatx0=1.0,x1;inti;
x1=（x0*x0*x0-1）/3;
printf（"\n%1.5f",x1）;
for（i=1;i<=50&&fabs（x1-x0）>1.0E-5;i++）
{
x0=x1;
x1=（x0*x0*x0-1）/3;
printf（"%1.5f",x1）;
}
printf（"\nyidaitimes:
%d",i）;
}
即当x0=2.0时其收敛，迭代次数为7次，结果为：
而当x0=2.0时的输出结果为：
说明当x0=2.0时发散，无法用此初值求的根。
3.x=（3*x+1）1/3c语言编程为
#include
main（）
{
floatx0=1.0,x1;inti;
x1=pow（（3*x0+1）,3）;
printf（"\n%1.5f",x1）;
for（i=1;i<=50&&fabs（x1-x0）>1.0E-5;i++）
{
x0=x1;
x1=pow（（3*x0+1）,3）;
printf（"%1.5f",x1）;
}
printf（"\nyidaitimes:
%d",i）;
}
当x0=1.0时发散，其结果为：
4.x=1/（x*x-3）,c语言编程为:
#include
main（）
{
floatx0=1.0,x1;inti;
x1=1/（x0*x0-3）;
printf（"\n%1.5f",x1）;
for（i=1;i<=50&&fabs（x1-x0）>1.0E-5;i++）
{
x0=x1;
x1=1/（x0*x0-3）;
printf（"%1.5f",x1）;
}
printf（"\nyidaitimes:
%d",i）;
}
即当x0=1.0时收敛，迭代次数为6次结果为:
5.x=（3+1/x）1/2，c语言编程为：
#include
main（）
{
floatx0=1.0,x1;inti;
x1=sqrt（3+1/x0）;
printf（"\n%1.5f",x1）;
for（i=1;i<=50&&fabs（x1-x0）>1.0E-5;i++）
{
x0=x1;
x1=sqrt（3+1/x0）;
printf（"%1.5f",x1）;
}
printf（"\nyidaitimes:
%d",i）;
}
即当x0=1.0时收敛，迭代次数为6次，结果为：
6.x=x-（1/3）*（x*x*x-3*x-1）/（x*x-1），c语言编程为
#include
main（）
{
floatx0=1.2,x1;inti;
x1=x0-（1/3）*（x0*x0*x0-3*x0-1）/（x0*x0-1）;
printf（"\n%1.5f",x1）;
for（i=1;i<=50&&fabs（x1-x0）>1.0E-5;i++）
{
x0=x1;
x1=x0-（1/3）*（x0*x0*x0-3*x0-1）/（x0*x0-1）;
printf（"%1.5f",x1）;
}
printf（"\nyidaitimes:
%d",i）;
}
当x0=1.2时结果为：
6式其实是一个恒等式不能用于迭代运算！
二．结果讨论
1.对于同一函数的不同迭代格式会有不同的敛散情况。
2.对于取不同的初值会影响到函数求根的迭代次数，而且会影响到迭代格式的敛散性（如弟2个函数）。
3.编程时可以用终止准则xk+1-xk的绝对值4.对于什么样的迭代格式符合要求，书上已经给出结论：
假设f（x）在区间[a,b]可导，且满足条件：
（1）a<=f（x）<=b,x在[a,b]内。
（2）存在正数L《1，使得对任意的在[a,b]内的x都有f（x）的导数的绝对值<=L<1,则方程x=f（x）在[a,b]上存在唯一根，且对任意在[a,b]上的初始值x0，由迭代格式x=f（x）产生的迭代序列均收敛于此根。
（证明见书）
5.通过此次上机实验，应该认识到，通过计算机编程可以用简单迭代法求出非线性方程的跟，其中迭代格式和初值的选取对能否求出根及能否快速求出根起着重要的作用。
升级会员 