数值分析实验报告.docx

数值分析实验报告.docx

《数值分析实验报告.docx》由会员分享，可在线阅读，更多相关《数值分析实验报告.docx（19页珍藏版）》请在冰豆网上搜索。

数值分析实验报告

实验报告一

题目：

非线性方程求解

摘要：

非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。

本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。

利用二分法求解给定非线性方程的根，在给定的范围内，假设f（x,y）在[a,b]上连续，f（a）xf（b）<0且f（x）在（a,b）内仅有一实根，取中点一步步求得根的近似解，在满足精度要求时，即为所求的值。

Newton法乃利用

，由递推产生近似于真值的解，但Newton法的初值选择好坏直接影响迭代的次数甚至迭代的收敛与发散。

即若x0偏离所求根较远，Newton法可能发散的结论。

并且本实验中还利用利用改进的Newton法求解同样的方程，且将结果与Newton法的结果比较分析。

前言：

（目的和意义）

掌握二分法与Newton法的基本原理和应用。

掌握二分法的原理，验证二分法，在选对有根区间的前提下，必是收敛，但精度不够。

熟悉Matlab语言编程，学习编程要点。

体会Newton使用时的优点，和局部收敛性，而在初值选取不当时，会发散。

数学原理：

对于一个非线性方程的数值解法很多。

在此介绍两种最常见的方法：

二分法和Newton法。

对于二分法，其数学实质就是说对于给定的待求解的方程f（x），其在[a,b]上连续，f（a）f（b）<0，且f（x）在[a,b]内仅有一个实根x*，取区间中点c，若，则c恰为其根，否则根据f（a）f（c）<0是否成立判断根在区间[a,c]和[c,b]中的哪一个，从而得出新区间，仍称为[a,b]。

重复运行计算，直至满足精度为止。

这就是二分法的计算思想。

Newton法通常预先要给出一个猜测初值x0，然后根据其迭代公式

产生逼近解x*的迭代数列{xk}，这就是Newton法的思想。

当x0接近x*时收敛很快，但是当x0选择不好时，可能会发散，因此初值的选取很重要。

另外，若将该迭代公式改进为

其中r为要求的方程的根的重数，这就是改进的Newton法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton法快的多。

程序设计：

本实验采用Matlab的M文件编写。

其中待求解的方程写成function的方式，如下

functiony=f（x）;

y=-x*x-sin（x）;

写成如上形式即可，下面给出主程序。

二分法源程序：

clear

%%%给定求解区间

b=1.5;

a=0;

%%%误差

R=1;

k=0;%迭代次数初值

while（R>5e-6）;

c=（a+b）/2;

iff12（a）*f12（c）>0;

a=c;

else

b=c;

end

R=b-a;%求出误差

k=k+1;

end

x=c%给出解

Newton法及改进的Newton法源程序：

clear

%%%%输入函数

f=input（'请输入需要求解函数>>','s'）

%%%求解f（x）的导数

df=diff（f）;

%%%改进常数或重根数

miu=2;

%%%初始值x0

x0=input（'inputinitialvaluex0>>'）;

k=0;%迭代次数

max=100;%最大迭代次数

R=eval（subs（f,'x0','x'））;%求解f（x0），以确定初值x0时否就是解

while（abs（R）>1e-8）

x1=x0-miu*eval（subs（f,'x0','x'））/eval（subs（df,'x0','x'））;

R=x1-x0;

x0=x1;

k=k+1;

if（eval（subs（f,'x0','x'））<1e-10）;

break

end

ifk>max;%如果迭代次数大于给定值，认为迭代不收敛，重新输入初值

ss=input（'mayberesultiserror,chooseanewx0,y/n?

>>','s'）;

ifstrcmp（ss,'y'）

x0=input（'inputinitialvaluex0>>'）;

k=0;

else

break

end

end

end

k;%给出迭代次数

x=x0；%给出解

结果分析和讨论：

1.用二分法计算方程

在[1，2]内的根。

（

下同）

计算结果为

x=1.40441513061523；

f（x）=-3.797205105904311e-007；

k=18；

由f（x）知结果满足要求，但迭代次数比较多，方法收敛速度比较慢。

2.用二分法计算方程

在[1，1.5]内的根。

计算结果为

x=1.32471847534180；

f（x）=2.209494846194815e-006；

k=17；

由f（x）知结果满足要求，但迭代次数还是比较多。

3.用Newton法求解下列方程

a）

x0=0.5；

计算结果为

x=0.56714329040978；

f（x）=2.220446049250313e-016；

k=4；

由f（x）知结果满足要求，而且又迭代次数只有4次看出收敛速度很快。

b）

x0=1；

c）

x0=0.45,x0=0.65；

当x0=0.45时，计算结果为

x=0.49999999999983；

f（x）=-8.362754932994584e-014；

k=4；

由f（x）知结果满足要求，而且又迭代次数只有4次看出收敛速度很快，实际上该方程确实有真解x=0.5。

当x0=0.65时，计算结果为

x=0.50000000000000；

f（x）=0；

k=9；

由f（x）知结果满足要求，实际上该方程确实有真解x=0.5，但迭代次数增多，实际上当取x0〉0.68时，x≈1，就变成了方程的另一个解，这说明Newton法收敛与初值很有关系，有的时候甚至可能不收敛。

4.用改进的Newton法求解，有2重根，取

x0=0.55；并与3.中的c）比较结果。

当x0=0.55时，程序死循环，无法计算，也就是说不收敛。

改

时，结果收敛为

x=0.50000087704286；

f（x）=4.385198907621127e-007；

k=16；

显然这个结果不是很好，而且也不是收敛至方程的2重根上。

当x0=0.85时，结果收敛为

x=1.00000000000489；

f（x）=2.394337647718737e-023；

k=4；

这次达到了预期的结果，这说明初值的选取很重要，直接关系到方法的收敛性，实际上直接用Newton法，在给定同样的条件和精度要求下，可得其迭代次数k=15，这说明改进后的Newton法法速度确实比较快。

结论：

对于二分法，只要能够保证在给定的区间内有根，使能够收敛的，当时收敛的速度和给定的区间有关，二且总体上来说速度比较慢。

Newton法，收敛速度要比二分法快，但是最终其收敛的结果与初值的选取有关，初值不同，收敛的结果也可能不一样，也就是结果可能不是预期需要的结果。

改进的Newton法求解重根问题时，如果初值不当，可能会不收敛，这一点非常重要，当然初值合适，相同情况下其速度要比Newton法快得多。

在编制程序过程中，需事先规定迭代的次数，若超过这个次数，还不收敛，则停止迭代，另选初值。

实验报告二

题目：

Romberg积分法

摘要：

应用数值方法进行求解积分问题已经有着很广泛的应用，本实验基于Romberg积分法来解决一类积分问题。

Romberg积分法的算法基础是基于复化梯形公式的，T0,K将区间进行2k等份的复化梯形求积公式的计算结果，根据理查逊外推加速方法推出Tm,K来。

很明显，Romberg积分法为高速有效，易于编制程序，非常适合于计算机计算的方法。

另外，Romberg积分法计算过程在计算机上实现时，只需开辟一个一维数组，即每次计算的结果Tm,K，可存放在T0,K位置上，其最终结果Tm,0存放在T0,0位置上。

前言：

（目的和意义）

1.理解和掌握Romberg积分法的原理；

2.学会使用Romberg积分法；

3.明确Romberg积分法的收敛速度及应用时容易出现的问题；

4．体会算法的思想，体会计算机在此类积分中所发挥出来的巨大的优势；

5．学会用matlab语言编制程序，解决比较繁琐的积分问题。

数学原理：

考虑积分

，欲求其近似值，通常有复化的梯形公式、Simpsion公式和Cotes公式。

但是给定一个精度，这些公式达到要求的速度很缓慢。

如何提高收敛速度，自然是人们极为关心的课题。

为此，记T1,k为将区间[a,b]进行2k等分的复化的梯形公式计算结果，记T2,k为将区间[a,b]进行2k等分的复化的Simpsion公式计算结果，记T3,k为将区间[a,b]进行2k等分的复化的Cotes公式计算结果。

根据Richardson外推加速方法，可以得到收敛速度较快的Romberg积分法。

其具体的计算公式为：

1.准备初值，计算

2.按梯形公式的递推关系，计算

3.按Romberg积分公式计算加速值

m=1,2,…，k=i-m

4.精度控制。

对给定的精度

，若

则终止计算，并取

为所求结果；否则返回2重复计算，直至满足要求的精度为止。

程序设计：

本实验采用Matlab的M文件编写。

其中待积分的函数写成function的方式，例如如下

functionyy=f（x,y）;

yy=x.^3;

写成如上形式即可，下面给出主程序

Romberg积分法源程序

%%%Romberg积分法

clear

%%%积分区间

b=3;

a=1;

%%%精度要求

R=1e-5;

%%%应用梯形公式准备初值

T（1,1）=（b-a）*（f（b）+f（a））/2;

T（1,2）=T（1,1）/2+（b-a）/2*f（（b+a）/2）;

T（2,1）=（4*T（1,2）-T（1,1））/（4-1）;

j=2;

m=2;

%%%主程序体%%%

while（abs（T（m,1）-T（m-1,1））>R）;%%%精度控制

j=j+1;

s=0;

forp=1:

2^（j-2）;

s=s+f（a+（2*p-1）*h/（2^（j-1）））;

end

T（1,j）=T（1,j-1）/2+h*s/（2^（j-1））;%%%梯形公式应用

form=2:

j;

k=（j-m+1）;

T（m,k）=（（4^（m-1））*T（m-1,k+1）-T（m-1,k））/（4^（m-1）-1）;

end

end

%%%给出Romberg积分法的函数表

I=T（m,1）

结果分析和讨论：

进行具体的积分时，精度取R=1e-8。

1.求积分

。

精确解I=24999676。

运行程序得Romberg积分法的函数表为

1.0e+007*

4.701015200000003.050229500000002.63753307500000

2.499967600000002.499967600000000

2.4999676000000000

由函数表知Romberg积分给出的结果为2.4999676*10^7，与精确没有误差，

精度很高。

2.求积分

。

精确解I=ln3=1.09861228866811。

运行程序得Romberg积分法的函数表为

1.333333333333331.166********6671.116666666666671.103210678210681.099767701563031.098901515168461.09868461878559

1.111111111111111.100000000000001.098725348725351.098620042680481.098612786370271.098612319991300

1.099259259259261.098640371973711.098613022277491.098612302616251.0986122888993700

1.098630548366001.098612588155331.098612291193061.09861228868164000

1.098612517723131.098612290028501.098612288671790000

1.098612289805931.0986122886704600000

1.09861228867019000000

从积分表中可以看出程序运行的结果为1.09861228867019，取8位有效数字，满足要求。

3.求积分

。

直接按前面方法进行积分，会发现系统报错，出现了0为除数的现象。

出现这种情况的原因就是当x=0时，被积函数分母出现了0，如果用一个适当的小数

（最好不要小于程序给定的最小误差值，但是不能小于机器的最大精度）来代替，可以避免这个问题。

本实验取

，可得函数表为：

0.920735483196590.939793275001900.944513511714170.945690853594890.94598501993743

0.946145872270340.946086923951600.946083300888460.946083075384950

0.946082994063680.946083059350920.9460830603513800

0.946083060387220.94608306036726000

0.946083060367180000

故该函数的积分为0.94608306036718，取8位有效数字。

4.求积分

本题的解析解很难给出，但运用Romberg积分可以很容易给出近似解，函数表为：

0.420735492403950.334069725829240.315975360759220.311680239480940.310620366809490.31035626065456

0.305181136971000.309943905735880.310248532388180.310267075919000.310268225269590

0.310261423653540.310268840831670.310268312154390.3102683018929600

0.310268958564650.310268303762690.31026830173008000

0.310268301194840.310268301722110000

0.3102683017226200000

故该函数的积分为0.31026830172262，取8位有效数字。

结论：

应用数值方法进行求解积分问题已经有着很广泛的应用，实验发现，Romberg积分通常要求被积函数在积分区间上没有奇点。

如有奇点，且奇点为第一间断点，那么采用例3的方法，还是能够求出来的，否则，必须采用其它的积分方法。

当然，Romberg积分的收敛速度还是比较快的，易于编程控制，适合于计算机计算。

但它也有一个主要缺点，每当把区间对分后，就要对被积函数f（x）计算它在新分点处的值，而这些函数值的个数也是成倍增加的。

另外，Romberg积分法计算过程在计算机上实现时，只需开辟一个一维数组，即每次计算的结果Tm,K，可存放在T0,K位置上，其最终结果Tm,0存放在T0,0位置上，节省不少内存。

实验报告三

题目：

常微分方程初值问题的数值解法

摘要：

常微分方程初值问题的数值解法就是求初值问题的解在点列xn=a+nh上的近似值yn，其中h称为步长。

常微分方程初值问题的数值解法一般分为两大类，即单步法和多步法。

本实验主要采用经典四阶的P-K方法和四阶Adams预测-校正方法来求解常微分方程的数值解。

发现P-K方法的优点是单步法、精度高，在计算过程中便于改变步长，缺点是计算量较大，每前进一步需要计算四次函数值f（x,y）。

与P-K方法相比，Adams可以节省计算量，每前进一步需要计算二次函数值f（x,y）；缺点是它不是自开始的，需要先知道前面四个点的值y0，y1，y2，y3，因此，它不能独立使用，另外，它也不便于改变步长。

前言：

（目的和意义）

通过编写程序，进行上机计算，使得对常微分方程初值问题的数值解法有更深刻的理解，掌握单步法和线性多步法是如何进行实际计算的及两类方法的适用范围和优缺点，特别是对这两类方法中最有代表性的方法:

P-K方法和Adams方法及预测-校正方法有更好的理解。

通过这两种方法的配合使用，掌握不同的方法如何配合在一起，解决实际问题。

数学原理：

对于一阶常微分方程初值问题

（1）

的数值解法是近似计算中很中的一部分。

常微分方程的数值解法通常就是给出定义域上的n个等距节点，求出所对应的函数值yn。

通常其数值解法可分为两大类：

1.单步法：

这类方法在计算yn+1的值时，只需要知道xn+1、xn和yn即可，就可算出。

这类方法典型有欧拉法和P-K法。

2.多步法：

这类方法在计算yn+1的值时，除了需要知道xn+1、xn和yn值外，还需要知道前k步的值。

典型的方法如Adams法。

经典的P-K法是一个四阶的方法。

它最大的优点就是它是单步法，精度高，计算过程便于改变步长。

其缺点也很明显，计算量大，每前进一步就要计算四次函数值f。

它的具体的计算公式如式

（2）所示。

四阶Adams预测-校正方法是一个线性多步法，它是由Adams显式公式

（2）

和隐式公式组成，其计算公式如式（3）所示。

预测

（3a）

求导

（3b）

校正

（3c）

求导

（3d）

将局部截断误差用预测值和校正值来表示，在预测和校正的公式中分别以它们各自的阶段误差来进行弥补，可期望的到精度更高的修正的预测-校正公式为：

预测

修正

求导

校正

修正

求导

由于开始时无预测值和校正值可以利用，故令p0=c0=0，以后按上面进行计算。

此方法的优点是可以减少计算量；缺点是它不是自开始的，需要先知道前面的四个点的值

，因此不能单独使用。

另外，它也不便于改变步长。

程序设计：

本实验采用Matlab的M文件编写。

其中待求的微分方程写成function的方式，如下

functionyy=g（x,y）;

yy=-x*x-y*y;

写成如上形式即可，下面给出主程序。

经典四阶的P-K方法源程序

clear

%%%步长选取

h=0.1;

%%%初始条件，即x=0时，y=1。

y

（1）=1;

%%%求解区间

a=0;

b=2;

%%%迭代公式

forx=a:

h:

b-h;

k1=g（x,y（（x-a）/h+1））;%%

，下同。

k2=g（x+h/2,y（（x-a）/h+1）+h/2*k1）;

k3=g（x+h/2,y（（x-a）/h+1）+h/2*k2）;

k4=g（x+h,y（（x-a）/h+1）+h*k3）;

y（（x-a）/h+2）=y（（x-a）/h+1）+h*（k1+2*k2+2*k3+k4）/6;

end

四阶Adams预测-校正方法源程序

%%%步长选取

h=0.1;

%%%初始条件

y

（1）=1;

%%%求解区间

a=0;

b=2;

%%%应用RK迭代公式计算初始值y0,y1,y2,y3

forx=a:

h:

a+2*h;

k1=g（x,y（（x-a）/h+1））;

k2=g（x+h/2,y（（x-a）/h+1）+h/2*k1）;

k3=g（x+h/2,y（（x-a）/h+1）+h/2*k2）;

k4=g（x+h,y（（x-a）/h+1）+h*k3）;

y（（x-a）/h+2）=y（（x-a）/h+1）+h*（k1+2*k2+2*k3+k4）/6;

end

%%%应用预测校正法求解

c（4）=0;%%%校正初值

p（4）=0;%%%预测初值

f

（1）=g（a+0*h,y

（1））;

，且将该值存在数组f中。

f

（2）=g（a+1*h,y

（2））;

f（3）=g（a+2*h,y（3））;

f（4）=g（a+3*h,y（4））;

forn=4:

（b-a）/h;

%%%P预测

p（n+1）=y（n）+h/24*（55*f（n）-59*f（n-1）+37*f（n-2）-9*f（n-3））;

%%%M修正

m（n+1）=p（n+1）+251/270*（c（n）-p（n））;

%%%E求导

f（n+1）=g（a+（n+1-1）*h,m（n+1））;

%%%C校正

c（n+1）=y（n）+h/24*（9*f（n+1）+19*f（n）-5*f（n-1）+f（n-2））;

%%%M修正

y（n+1）=c（n+1）-19/270*（c（n+1）-p（n+1））;

%%%E求导

f（n+1）=g（a+（n+1-1）*h,y（n+1））;

end

结果分析和讨论：

计算实例一：

对初值问题

取步长h=0.1;计算在[-1，0]上的数值解。

在计算机上，运用所编的程序进行了运算，结果如下，其中第一列为在区间上的等分点，第二列为运用P-K法的计算结果，第三列为Adams预测-校正法计算结果。

由于本题的解析解很难求出，无法看出精度如何，为此进行第二实例计算。

第一实例计算结果

-1.0000000000000000

-0.900000000000000.090047357462400.09004735746240

-0.800000000000000.160726883901280.16072688390128

-0.700000000000000.213482650382680.21348265038268

-0.600000000000000.250368369544590.25036768670937

-0.500000000000000.273775247604510.27377644362732

-0.400000000000000.286221918318140.28622489850956

-0.300000000000000.290213341573770.29021772444822

-0.200000000000000.288160170936890.28816545338237

-0.100000000000000.282343663961140.28234937949911

00.274910519234620.27491630159737

计算实例二：

对初值问题

取步长h=0.1;计算在[0，1.5]上的数值解。

本题的解析解为

。

同样在计算机上，运用所编的程序进行了运算，结果如下，其中第一列为在区间上的等分点，第二列为运用P-K法的计算结果，第三列为Adams预测-校正法计算结果，第四列为精确解。

第二实例计算结果

01.000000000000001.000000000000001.00000000000000

0.100000000000001.095445531693091.095445531693091.09544511501