数值计算大作业——刘峰Word文档格式.doc

数值计算大作业——刘峰Word文档格式.doc

《数值计算大作业——刘峰Word文档格式.doc》由会员分享，可在线阅读，更多相关《数值计算大作业——刘峰Word文档格式.doc（12页珍藏版）》请在冰豆网上搜索。

采用牛顿迭代法，牛顿迭代法的数学原理是，对于方程，如果是线性函数，则它的求根是很容易的，牛顿迭代法实质上是一种线性化方法，其基本思想是将非线性方程逐步归结为某种线性方程来求解。

设已知方程有近似根（假定），将函数在点进行泰勒展开，有

于是方程可近似地表示为

这是个线性方程，记其根为，则的计算公式为

，

这就是牛顿迭代法，简称牛顿法。

3.程序设计

作出函数的图像，大概估计出根的位置

fplot（'

1000*exp（x）+（435*x）*（exp（x）-1）-1564'

[03]）;

grid

大概估计出初始值x=0.5

function[p1,err,k,y]=newton（f,df,p0,delta,max1）

%f是非线性系数

%df是f的微商

%p0是初始值

%dalta是给定允许误差

%max1是迭代的最大次数

%p1是牛顿法求得的方程近似解

%err是p0误差估计

%k是迭代次数

p0,feval（'

f'

p0）

fork=1:

max1

p1=p0-feval（'

p0）/feval（'

df'

p0）;

err=abs（p1-p0）;

p0=p1;

p1,err,k,y=feval（'

p1）

if（err<

delta）|（y==0）,

break,end

p1,err,k,y=feval（'

end

functiony=f（x）

y=1000000*exp（x）+435000*（exp（x）-1）/x-1564000;

functiony=df（x）

y=1000000*exp（x）+435000*（exp（x）/x-（exp（x）-1）/x^2）;

4.结果分析与讨论

在MATLAB中的commandwindow输入

newton（'

'

1.2,10^（-4）,10）

运行后得出结果

p0=0.5000

p1=0.1679err=0.3321k=1y=9.2415e+004

p1=0.1031err=0.0648k=2y=2.7701e+003

p1=0.1010err=0.0021k=3y=2.6953

p1=0.1010err=2.0129e-006k=4y=2.5576e-006

ans=0.1010

运算后的结果为，通过这个数值来预测第二年年末的人口数，

t=2时候对于

实践表明，当初始值难以确定时，迭代法就不一定收敛了，因此要根据问题实际背景或者二分法先得一个较好的初始值，然后再进行迭代；

再者迭代函数选择不合适的话，采用不动点迭代法也有可能出现不收敛的情况；

因此我采用的是牛顿法。

题目二：

线性方程组求解

假设一个物体可以位于个等距点的任意位置，当物体在位置时，它只能等可能的移动到或者，而不能直接移动到其他任何位置，概率表示物体从位置开始在到达右端点之前到达左端点的概率，显然，且有

既有下面方程组：

取对方程组进行求解（迭代法或者直接法）。

在解微分方程的边值问题、热传导方程以及船体数学放样中建立的三次样条函数等工程技术问题时，经常遇到下面形式的线性方程组：

=

方程简记，该线性方程称为三对角线方程组，其系数矩阵A满足条件

所以为弱对角阵可以采用追赶法进行计算，利用三对角矩阵的LU分解建立计算量更少的线性方程组求解公式。

将系数矩阵A进行克劳特分解，即A分解为下三角矩阵和单位上三角矩阵的乘积；

A==

其中，，为待定系数，直接利用矩阵乘法公式可得

，，

于是推得计算，，的公式

，；

，，；

由此计算出L和U中的全部元素，完成了系数矩阵A的克劳特分解。

求解线性方程组等价于求解和。

因而得到解三对角线性方程组的追赶法公式

（1）计算的递推公式：

（2）解

（3）解

我们将计算系数和称为追的过程，将计算方程组的解称为赶的过程。

整个过程为追赶法的思想。

functionx=chase（a,b,c,f）

%求解线性方程组Ax=f，其中A是三对角阵

%a是矩阵A的下对角线元素a

（1）=0

%b是矩阵A的对角线元素

%c是矩阵A的上对角线元素c（N）=0

%f是方程组的右端向量

n=length（b）;

ifn-1==length（a）

fori=n-1:

-1:

1

a（i+1）=a（i）;

end

end

c

（1）=c

（1）/b

（1）;

f

（1）=f

（1）/b

（1）;

fori=2:

n-1

b（i）=b（i）-a（i）*c（i-1）;

c（i）=c（i）/b（i）;

f（i）=（f（i）-a（i）*f（i-1））/b（i）;

f（n）=（f（n）-a（n）*f（n-1））/（b（n）-a（n）*c（n-1））;

fori=n-1:

1f（i）=f（i）-c（i）*f（i+1）;

x=f;

A的系数矩阵为

A=[1,-0.5,0,0,0,0,0,0,0,0;

-0.5,1,-0.5,0,0,0,0,0,0,0;

0,-0.5,1,-0.5,0,0,0,0,0,0;

0,0,-0.5,1,-0.5,0,0,0,0,0;

...

0,0,0,-0.5,1,-0.5,0,0,0,0;

0,0,0,0,-0.5,1,-0.5,0,0,0;

0,0,0,0,0,-0.5,1,-0.5,0,0;

0,0,0,0,0,0,-0.5,1,-0.5,0;

0,0,0,0,0,0,0,-0.5,1,-0.5;

0,0,0,0,0,0,0,0,-0.5,1;

]

所以在MATLAB命令窗口输入

>

a=[-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,0]

b=[1,1,1,1,1,1,1,1,1,1]

c=[-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,0]

f=[0.5,0,0,0,0,0,0,0,0,0]

得到此题中的a,b,c,f矩阵：

a=

-0.5000-0.5000-0.5000-0.5000-0.5000-0.5000-0.5000-0.5000-0.5000

b=

1111111111

c=

-0.5000-0.5000-0.5000-0.5000-0.5000-0.5000

-0.5000-0.5000-0.50000

f=

0.5000000000000

然后在MATLAB中调用之前保存的迭代法函数function，在命令窗口中输入：

chase（a,b,c,f）

回车得到结果：

x=chase（a,b,c,f）

x=

0.90000.80000.70000.60000.50000.40000.30000.20000.10000

追赶法为一种特殊的LU分解法。

追赶法是求解三对角矩阵的常用方法，但从整体编程角度分析，其程序编写较迭代法复杂，但通用性较好。

追赶法求解三对角矩阵不但节省存储单元，而且可以减少计算量，是工程技术中比较常用的数学工具。

三、数值积分

1、题目

卫星轨道是一个椭圆，椭圆周长的计算公式是,这里是椭圆的半长轴,是地球中心与轨道中心（椭圆中心）的距离,记为近地点距离,为远地点距离,公里为地球半径,则,某人造卫星近地点距离公里,远地点距离公里,试用Romberg方法求卫星轨道的周长，精确到。

龙贝格方法是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。

作为一种外推算法,它在不增加计算量的前提下提高了误差的精度。

龙贝格方法的主要过程是将粗糙的梯形公式逐步加工成精度较高的辛普森公式和科特斯公式的方法称为龙贝格方法。

复化梯形公式

在复化梯形公式中，每个内节点既是前一个小区间的终点，又是后一个小区间的起点，因此上式可以改写为

复化梯形公式余项

复化梯形公式的递推公式为

复化辛普森求积公式

与复化梯形公式类似，每个内节点需用两次，因此有

显然复化辛普森公式在n趋于无穷大时，他的收敛速度比复化梯形公式更快。

以表示二分k次后求得的梯形值，且以表示序列的m次加速度，理查森外推法的递推公式可写成

龙贝格算法的计算过程如下：

（1）取求

（2）利用变步长梯形公式，其中k为区间的二分次数，即

或

（3）依横行次序求加速值，逐个求出的第k行其余各元素

（4）当相邻对角元素之差的绝对值小于预先给定的精度时，终止计算。

表3-1龙贝格算法递推表

k

h

b-a

1

2

3

4

functionR=romberg（f,a,b,n）

formatlong

R=zeros（[n+1,n+1]）;

R（0+1,0+1）=（b-a）/2*（feval（f,a）+feval（f,b））;

fori=1:

n,h=（b-a）/2^i;

s=0;

fork=1:

2^（i-1）,

s=s+feval（f,a+（2*k-1）*h）;

end

R（i+1,0+1）=R（i-1+1,0+1）/2+h*s