非线性方程的数值计算方法实验解析.docx
《非线性方程的数值计算方法实验解析.docx》由会员分享，可在线阅读，更多相关《非线性方程的数值计算方法实验解析.docx（14页珍藏版）》请在冰豆网上搜索。
非线性方程的数值计算方法实验解析
非线性方程的数值计算方法实验
一、实验描述：
在科学研究和工程实践中，经常需要求解大量的非线性方程。
本实验正是通过计算机的程序设计，使用迭代法、波尔查诺二分法、试值法、牛顿-拉夫森法和割线法，来实现非线性方程的求解。
本实验中通过对各种方法的实践运用，可以比较出各种方法的优缺点。
并且，通过完成实验，可加深对各种方法的原理的理解，熟悉掌握C语言在这些方法中的运用。
二、实验内容：
1、求函数
的不动点（尽可能多）近似值，答案精确到小数点后12位；
2、如果在240个月内每月付款300美元，求解满足全部年金A为500000美元的利率I，的近似值（精确到小数点后10位）。
3、利用加速牛顿-拉夫森算法，用其求下列函数M阶根p的近似值。
（a）、f（x）=（x-2）5，M=5,p=2，初始值p0=1。
（b）、f（x）=sin（x3），M=3,p=0，初始值p0=1。
（c）、f（x）=（x-1）ln（x）,M=2,p=1，初始值p0=2。
4、设投射体的运动方程为：
y=f（t）=9600（1-e-t/15）-480t
x=r（t）=2400（1-e-t/15）
（a）求当撞击地面时经过的时间，精确到小数点后10位。
（b）求水平飞行行程，精确到小数点后10位。
三、实验原理：
（1）、不动点迭代法：
它是一种逐次逼近的方法,即用某个固定公式反复校正根的近似值,使之逐步精确化,最后得到满足精度要求的结果。
它利用计算机运算速度快，适合做重复性操作的特点，让计算机对一个函数进行重复执行，在每次执行这个函数时，都从变量的原值推出它的一个新值，直至推出最终答案为止。
迭代法一般可用于寻找不动点，即：
存在一个实数P，满足P=g（P），则称P为函数g（x）的一个不动点。
且有定理：
若g（x）是一个连续函数，且
是由不动点迭代生成的序列。
如果
=P,则P是g（x）的不动点。
所以，不动点的寻找多用迭代法。
（2）、波尔查诺二分法：
起始区间[a,b]必须满足f（a）与f（b）的符号相反的条件。
由于连续函数y=f（x）的图形无间断，所以它会在零点x=r处跨过x轴，且r在区间内。
通过二分法可将区间内的端点逐步逼近零点，直到得到一个任意小的包含零点的间隔。
二分法定理：
设f
C（a,b）,且存在数r
[a,b]满足f（r）=0。
如果f（a）和f（b）的符号相反，且
为二分法生成的中点序列，则：
其中n=0,1,…
（1）
这样，序列
收敛到零点x=r即可表示为：
（2）
（3）、试值法：
假设一个函数中，有f（a）和f（b）符号相反。
二分法使用区间[a,b]的中点进行下一次迭代。
如果找到经过点（a,f（a））和（b,f（b））的割线L与x轴的交点（c,0）,则可得到一个更好的近似值。
为了寻找值c，定义了线L的斜率m的两种表示方法，一种表示方法为：
m=
（3）
这里使用了点（a,f（a））和（b,f（b））。
另一种表示方法为：
m=
（4）
这里使用了点（c,0）和（b,f（b））。
使式（3）和式（4）的斜率相等，则有：
（5）
为了更容易求解c，可进一步表示为：
c=b-
（6）
这样会出现3种可能性：
如果f（a）和f（c）的符号相反，则在[a,c]内有一个零点。
如果f（c）和f（b）的符号相反，则在[c,b]内有一个零点。
如果f（c）=0，则c是零点。
然后，可按二分法的方法进行下一步运算。
（4）、牛顿-拉夫森法：
此法根据，牛顿-拉夫森定理：
设f
[a,b],且存在数p
[a,b],
满足f（p）=0。
如果
（p）≠0,则存在一个数
>0，对任意初始近似值p0
[p-
p+
]，使得由如下迭代定义的序列
收敛到p:
=g（pk-1）=pk-1－
其中k=1,2,…（7）
其中，函数g（x）由如下定义：
g（x）=x--
（8）
且被称为牛顿-拉夫森迭代函数。
由于f（p）=0,显然g（p）=p。
这样，通过寻找函数的不动点，可以实现寻找方程f（x）=0的根的牛顿-拉夫森迭代。
附：
定理（牛顿-拉夫森迭代的加速收敛）：
设牛顿-拉夫森算法产生的序列线性收敛到M阶根x=p,其中M>1,则牛顿-拉夫森迭代公式：
=
-M
（9）
（5）、割线法：
割线法包含的公式与试值法的公式一样，只是在关于如何定义每个后续项的逻辑判定上不一样。
需要两个靠近点（p,0）的初始点（
f（
））和（
f（
））。
可根据两点迭代法公式，得到一般项:
=g（
）=
-
（10）
四、结果计算及分析：
1、函数g（x）=xx-cos（x）的不动点的迭代（迭代法）：
（a）计算结果：
（b）结果分析:
此题经过matlab图形仿真，可以看出其解的大致范围在[0.8,1.2]之间。
此结果是在取p0=0.879的情况下得出的，其误差精度取为0.000000000001。
从迭代过程的误差收敛速度可以看出不动点迭代的误差收敛较慢，所需迭代次数较多。
2、满足条件的利率I的计算（试值法）：
（a）、计算结果：
（b）、结果分析：
此题经过一定的计算分析，将书中所给公式结合题中的条件可得出最终的计算式：
[
-1]-500000=0,求解其中的I即为利率。
又经过一定的计算分析，取其初始计算区间为[0.1,0.2]，误差精度取为0.000000000001。
从其误差收敛的速度可以看出，试值法的误差收敛速度快，所需的迭代次数少。
3、利用加速牛顿-拉夫森算法，求函数f（x）=（x-2）5、f（x）=sin（x3）和f（x）=（x-1）ln（x）的M阶根p的近似值。
其中每个函数的M、p0和最终结果p已给出。
（a）、计算结果：
（b）、结果分析：
由于C语言函数库中没有求导功能函数，而此题所要求用的加速牛顿-拉夫森算法的计算公式
=
-M
中要用到题中所给函数的导函数，所以求出它们的导函数分别为
=5（x-2）2,
=3x2cos（x3）和
=In（x）-1-
。
此题的误差精度取为0.00000000001。
从输出结果中可以大致判断出，加速牛顿-拉夫森法的收敛速度非常快。
4、给出投射物体的运动方程y=f（t）=9600（1-e-t/15）-480t和x=r（t）=2400（1-e-t/15）求物体落地时经过的时间和水平飞行程。
（割线法）
（a）计算结果：
（b）、结果分析：
此题用了割线法的思路进行程序设计，由于割线法的公式要求，需要给出t0,t1两个初始值。
经过一定的分析，给出t0=5.0，t1=15.0，误差精度设为0.0000000001。
从输出结果中可以看出，割线法的误差收敛速度比较快，所需的迭代次数比较少。
五、实验结果分析：
经过以上四个实验的比较分析，结合书上的实例与原理分析，可以得出。
不动点迭代法比起其他的迭代法，迭代速度慢、次数多。
试值法的迭代次相对少一些，但若对其初始区间估算得太大的话，会导致迭代步数增多，甚至会导致发散。
割线法比试值法要快得多，可以达到相对较高的精度，且在其迭代过程中每步只需一次新的函数赋值，其收敛阶一般能达到1.1618。
加速牛顿-拉夫森法再求多重根时速度很快，收敛阶能达到2。
附件：
第一题：
#include
#include
#include//为了调用FLT_EPSILON，防止出现分母为零的情况。
//
#defineN1000//保证有足够多的运算次数能达到比较精确的结果，且便于修改。
//
#definepre0.0000000000001//使循环能够在达到某一误差精度后停下来。
//
doublemain（）
{
doublep[N],m,err[N],relerr[N];
intk,n;
p[0]=0.879;//给定一个合理的初始值以此来进行迭代。
//
for（k=0;k{
m=p[k]-cos（p[k]）;
p[k+1]=pow（p[k],m）;//进行迭代的方程，由书上给出。
//
err[k]=fabs（p[k+1]-p[k]）;//对绝对误差的运算。
//
relerr[k]=err[k]/（fabs（p[k+1]）+FLT_EPSILON）;//对相对误差的运算。
//
n=k;//把k赋给n,使下面输出时能够在合适的时候停下来。
//
if（err[k]//
}
for（k=0;kprintf（"\nX[%d]=%14.12ferr[%d]=%14.12frelerr[%d]=%14.12f",k,p[k],k,err[k],k,relerr[k]）;//对迭代点及误差的输出。
//
getch（）;
return0;
}
第二题：
#include
#include
#include
#defineN10000
#definepre0.000000000001//使循环能够在达到某一误差精度后停下来。
//
doublemain（）
{
doublef（doublex）;
doublea[N],b[N],I[N],err[N],relerr[N];
inti,n;
a[0]=0.10,b[0]=0.20;//预估它们的初始区间//
for（i=0;i{
I[i]=b[i]-（f（b[i]）*（b[i]-a[i]））/（f（b[i]）-f（a[i]））;//试值法公式，有书上给出。
//
if（（f（a[i]）*f（I[i]））>0）
a[i+1]=I[i],b[i+1]=b[i];//根据试值法，判断区间。
//
else
a[i+1]=a[i],b[i+1]=I[i];//根据试值法，判断区间。
//
err[i]=fabs（I[i]-I[i-1]）;//对绝对误差的运算。
//
relerr[i]=err[i]/（fabs（I[i]）+FLT_EPSILON）;//对相对误差的运算。
//
n=i;//将i值赋给n,使输出能够在合适的时候停下来。
//
if（err[i]//
}
err[0]=I[0]-0;//定义初始误差，便于下面输出。
//
relerr[0]=err[0]/（fabs（I[0]）+FLT_EPSILON）;//定义初始误差，便于下面输出。
//
for（i=0;iprintf（"\nI[%d]=%12.10ferr[%d]=%12.10frelerr[%d]=%12.10f",i,I[i],i,err[i],i,relerr[i]）;//输出利率的计算过程。
//
getch（）;
return0;
}
doublef（doublex）
{
return（（3600/x）*（pow（1+x/12,240）-1）-500000）;//此题的计算公式，由书上给出。
//
}
第三题：
#include
#include
#include//为了调用FLT_EPSILON，防止出现分母为零的情况。
//
#definepre0.00000000001//使循环能够在达到某一误差精度后停下来。
//
#defineN1000//迭代次数，以保证能得到最精确的值。
//
doublemain（）
{
doublefa（doublex）;//声明fa函数。
//
doublefb（doublex）;//声明fb函数。
//
doublefc（doublex）;//声明fc函数。
//
doublea[N],b[N],c[N],err[N];//定义结果数组，与误差数组。
//
intk,i,j,m;
a[0]=1.0,b[0]=1.0,c[0]=2.0;//赋初始值，由书上给出。
//
for（k=0;k{
a[k+1]=a[k]-5.0*fa（a[k]）;//加速牛顿-拉夫森公式，见书上P64。
下同。
//
err[k]=fabs（2.0-a[k+2]）;//对绝对误差的运算。
k加2，以保证下面输出时其有足够多的值输出。
//
i=k;//将循环停止时的k值赋给i，便于在下面输出时有恰当的值使输出停止。
下同。
//
if（err[k]}
for（k=0;k{
b[k+1]=b[k]-3.0*fb（b[k]）;
err[k]=fabs（0.0-b[k]）;//对绝对误差的运算。
//
j=k;
if（err[k]}
for（k=0;k{
c[k+1]=c[k]-2.0*fc（c[k]）;
err[k]=fabs（1.0-c[k]）;//对绝对误差的运算。
//
m=k;
if（err[k]}
if（i>j）j=i;if（m>j）m=j;//对i,m,j进行比较，并求出最大值，并用于输出。
//
for（k=0;kprintf（"\na[%d]=%12.10fb[%d]=%12.10fc[%d]=%12.10f",k,a[k],k,b[k],k,c[k]）;
getch（）;
return0;
}
doublefa（doublex）
{
return（（1.0/5.0）*（x-2.0））;//书上所给的计算要求的公式。
下同。
//
}
doublefb（doublex）
{
return（（1.0/（3.0*x*x））*tan（x*x*x））;
}
doublefc（doublex）
{
return（（（x-1.0）*log（x））/（log（x）+1.0-1.0/x））;
}
第四题：
#include
#include//方便各种数学函数的调用。
//
#include//为了调用FLT_EPSILON，防止出现分母为零的情况。
//
#defineN1000//保证有足够多的运算次数能达到比较精确的结果，且便于修改。
//
#definepre0.0000000001//使循环能够在达到某一误差精度后停下来。
//
doublemain（）
{
doublef（doublet）;//声明f函数。
//
doubler（doublet）;//声明r函数。
//
doublet[N],err[N],relerr[N],dif[N];//使结果能够以数组形式输出，便于分析。
//
charch1[]="t",ch2[]="dif",ch3[]="err",ch4[]="relerr";//对输出进行标注。
//
intk,i;//定义i，以便在达到所需精度之后，对所有结果以及对最终结果进行输出。
//
t[0]=5.0,t[1]=15.0;//给出任意合法的初始值。
//
err[0]=fabs（t[1]-t[0]）,relerr[0]=err[0]/fabs（t[1]+FLT_EPSILON）,dif[0]=t[1]-t[0];//对各个结果的初始值进行运算。
//
for（k=1;k{
t[k+1]=t[k]-（f（t[k]）*（t[k]-t[k-1]））/（f（t[k]）-f（t[k-1]））;//用割线法给出的公式进行运算，以求出最终结果。
//
dif[k]=t[k+1]-t[k];
err[k]=fabs（dif[k]）;//对绝对误差的运算。
//
relerr[k]=err[k]/（fabs（t[k+1]）+FLT_EPSILON）;//对相对误差的运算。
//
i=k+1;//在达到所需精度后，在往后加一位得到更精确的结果。
//
if（err[k]//
}
printf（"%s%s%s%s",ch1,ch2,ch3,ch4）;//对输出进行标记。
//
for（k=0;kprintf（"\n%15.10f%15.10f%15.10f%15.10f",t[k],dif[k],err[k],relerr[k]）;//输出各结果。
//
printf（"\n"）;//使输出间空出一行，便于观察。
//
printf（"\n当物体撞击地面时经过了%12.10f秒。
",t[i]）;//输出最终所需时间。
//
printf（"\n物体水平飞行了%15.10f米。
",r（t[i]））;//调用r函数，输出最终在水平方向的飞行行程。
//
getch（）;
return0;
}
doublef（doublet）
{
return（9600*（1-exp（-t/15））-480*t）;//物体在Y方向的路程与时间的关系，由书上给出。
//
}
doubler（doublet）
{
return（2400*（1-exp（-t/15）））;//物体在x方向的路程与时间的关系，由书上给出。
//
}
升级会员 