数值计算大作业——刘峰

1、课 程 设 计课程名称：数值分析设计题目：数值计算大作业学 号：S315070064姓 名：刘峰完成时间：2015年10月25日题目一、非线性方程求根1.题目 假设人口随时间和当时人口数目成比例连续增长，在此假设下人口在短期内的增长建立数学模型。（1）如果令表示在时刻的人口数目，表示固定的人口出生率，则人口数目满足微分方程，此方程的解为；（2）如果允许移民移入且速率为恒定的，则微分方程变成，此方程的解为；假设某地区初始有1000000人，在第一年有435000人移入，又假设在第一年年底该地区人口数量1564000人，试通过下面的方程确定人口出生率，精确到；且通过这个数值来预测第二年年末的人口数

2、，假设移民速度保持不变。2.数学原理采用牛顿迭代法，牛顿迭代法的数学原理是，对于方程，如果是线性函数，则它的求根是很容易的，牛顿迭代法实质上是一种线性化方法，其基本思想是将非线性方程逐步归结为某种线性方程来求解。设已知方程有近似根(假定)，将函数在点进行泰勒展开，有于是方程可近似地表示为这是个线性方程，记其根为，则的计算公式为，这就是牛顿迭代法，简称牛顿法。3.程序设计作出函数的图像，大概估计出根的位置fplot('1000*exp(x)+(435*x)*(exp(x)-1)-1564',0 3);grid大概估计出初始值x=0.5function p1,err,k,y=new

3、ton(f,df,p0,delta,max1)% f是非线性系数% df是f的微商% p0是初始值% dalta是给定允许误差% max1是迭代的最大次数% p1是牛顿法求得的方程近似解% err是p0误差估计% k是迭代次数p0,feval('f',p0)for k=1:max1 p1=p0-feval('f',p0)/feval('df',p0); err=abs(p1-p0); p0=p1; p1,err,k,y=feval('f',p1) if(err<delta)|(y=0), break,end p1,err,k

4、,y=feval('f',p1)endfunction y=f(x)y=1000000*exp(x)+435000*(exp(x)-1)/x-1564000;function y=df(x)y=1000000*exp(x)+435000*(exp(x)/x-(exp(x)-1)/x2);4.结果分析与讨论在MATLAB中的command window输入newton('f','df',1.2,10(-4),10)运行后得出结果p0 =0.5000 p1 =0.1679 err =0.3321 k =1 y =9.2415e+004p1 =0.10

5、31 err =0.0648 k =2 y =2.7701e+003p1 =0.1010 err =0.0021 k =3 y =2.6953p1 =0.1010 err =2.0129e-006 k =4 y = 2.5576e-006ans =0.1010运算后的结果为，通过这个数值来预测第二年年末的人口数，t=2时候对于实践表明，当初始值难以确定时，迭代法就不一定收敛了，因此要根据问题实际背景或者二分法先得一个较好的初始值，然后再进行迭代；再者迭代函数选择不合适的话，采用不动点迭代法也有可能出现不收敛的情况；因此我采用的是牛顿法。题目二：线性方程组求解1.题目 假设一个物体可以位于个等距

6、点的任意位置，当物体在位置时，它只能等可能的移动到或者，而不能直接移动到其他任何位置，概率表示物体从位置开始在到达右端点之前到达左端点的概率，显然，且有既有下面方程组：取对方程组进行求解（迭代法或者直接法）。2.数学原理 在解微分方程的边值问题、热传导方程以及船体数学放样中建立的三次样条函数等工程技术问题时，经常遇到下面形式的线性方程组： =方程简记，该线性方程称为三对角线方程组，其系数矩阵A满足条件所以为弱对角阵可以采用追赶法进行计算，利用三对角矩阵的LU分解建立计算量更少的线性方程组求解公式。将系数矩阵A进行克劳特分解，即A分解为下三角矩阵和单位上三角矩阵的乘积；A= 其中，为待定系数，直

7、接利用矩阵乘法公式可得，于是推得计算，的公式，；，；，；由此计算出L和U中的全部元素，完成了系数矩阵A的克劳特分解。求解线性方程组等价于求解和。因而得到解三对角线性方程组的追赶法公式(1) 计算的递推公式：(2) 解(3) 解我们将计算系数和称为追的过程，将计算方程组的解称为赶的过程。整个过程为追赶法的思想。3.程序设计function x=chase (a,b,c,f) %求解线性方程组Ax=f，其中A是三对角阵 %a 是矩阵A的下对角线元素 a（1）=0 %b 是矩阵A的对角线元素 %c 是矩阵A的上对角线元素 c（N）=0 %f 是方程组的右端向量 n=length(b); if n-1

8、=length(a) for i=n-1:-1:1 a(i+1)=a(i); end endc(1)=c(1)/b(1); f(1)=f(1)/b(1);for i=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);endf(n)=(f(n)-a(n)*f(n-1)/(b(n)-a(n)*c(n-1);for i=n-1:-1:1 f(i)=f(i)-c(i)*f(i+1); endx=f;4.结果分析与讨论A的系数矩阵为A=1,-0.5,0,0,0,0,0,0,0,0;-0.5,1,-0.5,0,

9、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>>

10、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 = 1 1 1 1 1 1 1 1 1 1c = -0.5000 -0.5000 -0.5000 -0.5000 -0.5000 -0.5000 -0.5000 -0.5000 -0.5000

11、 0f = 0.5000 0 0 0 0 0 0 0 0 0然后在MATLAB中调用之前保存的迭代法函数function，在命令窗口中输入：chase(a,b,c,f)回车得到结果：>> x=chase(a,b,c,f)x =0.9000 0.8000 0.7000 0.6000 0.5000 0.4000 0.3000 0.2000 0.1000 0追赶法为一种特殊的LU分解法。追赶法是求解三对角矩阵的常用方法，但从整体编程角度分析，其程序编写较迭代法复杂，但通用性较好。追赶法求解三对角矩阵不但节省存储单元，而且可以减少计算量，是工程技术中比较常用的数学工具。三、数值积分1、题目

12、卫星轨道是一个椭圆，椭圆周长的计算公式是, 这里是椭圆的半长轴, 是地球中心与轨道中心(椭圆中心)的距离, 记为近地点距离, 为远地点距离, 公里为地球半径,则, 某人造卫星近地点距离公里,远地点距离公里, 试用Romberg方法求卫星轨道的周长，精确到。2.数学原理龙贝格方法是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。 作为一种外推算法, 它在不增加计算量的前提下提高了误差的精度。 龙贝格方法的主要过程是将粗糙的梯形公式逐步加工成精度较高的辛普森公式和科特斯公式的方法称为龙贝格方法。复化梯形公式 在复化梯形公式中，每个内节点既是前一个小区间的终点，又

13、是后一个小区间的起点，因此上式可以改写为复化梯形公式余项 复化梯形公式的递推公式为复化辛普森求积公式与复化梯形公式类似，每个内节点需用两次，因此有显然复化辛普森公式在n趋于无穷大时，他的收敛速度比复化梯形公式更快。以表示二分k次后求得的梯形值，且以表示序列的m次加速度，理查森外推法的递推公式可写成龙贝格算法的计算过程如下：(1) 取求(2) 利用变步长梯形公式，其中k为区间的二分次数，即或(3) 依横行次序求加速值，逐个求出的第k行其余各元素(4) 当相邻对角元素之差的绝对值小于预先给定的精度时，终止计算。表3-1龙贝格算法递推表kh0b-a12343.程序设计function R=rombe

14、rg(f,a,b,n)format longR=zeros(n+1,n+1);R(0+1,0+1)=(b-a)/2*(feval(f,a)+feval(f,b);for i=1:n,h=(b-a)/2i; s=0; for k=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;endfor j=1:n,fac=1/(4j-1); for m=j:n, R(m+1,j+1)=R(m+1,j-1+1)+fac*(R(m+1,j-1+1)-R(m-1+1,j-1+1); endend4.结果分析与讨论本题根据算法原理在matlab中编写完龙贝格算法的自定义程序后，直接输入符合格式的函数积分就可得到相应轨道周长。调用MATLAB龙贝格算法的函数后可算得 R = romber