1、 课 程 设 计课程名称： 数值分析设计题目： 数值计算大作业学 号： S315070064姓 名： 刘峰完成时间： 2015 年 10 月 25 日- 1 -题目一、非线性方程求根1.题目 假设人口随时间和当时人口数目成比例连续增长，在此假设下人口在短期内的增长建立数学模型。（1）如果令 表示在 时刻的人口数目， 表示固定的人口出生率，则人口()Ntt数目满足微分方程 ，此方程的解为 ；()dt0()=tNte（2）如果允许移民移入且速率为恒定的 ，则微分方程变成 ，v()()dtNtv此方程的解为 ；0()=+(1)ttvNte假设某地区初始有 1000000 人，在第一年有 435000
2、 人移入，又假设在第一年年底该地区人口数量 1564000 人，试通过下面的方程确定人口出生率 ，精确到 ；且通过这个数值来预测第二年年末的人口数，假设移民速度 保持不410 v变。 435015640=(1)ee2.数学原理采用牛顿迭代法，牛顿迭代法的数学原理是，对于方程 ，如果0)(xf是线性函数，则它的求根是很容易的，牛顿迭代法实质上是一种线性化方)(xf法，其基本思想是将非线性方程 逐步归结为某种线性方程来求解。0)(xf设已知方程 有近似根 (假定 )，将函数 在点 进行0)(xfkf )(xfk泰勒展开，有 .)()(kkkxfxf于是方程 可近似地表示为0)(xf 0)()(kk
3、f这是个线性方程，记其根为 1kx，则 1的计算公式为- 2 -，)(1kkkxfx ,210这就是牛顿迭代法，简称牛顿法。3.程序设计作出函数的图像，大概估计出根的位置fplot(1000*exp(x)+(435*x)*(exp(x)-1)-1564,0 3);grid大概估计出初始值 x=0.5function p1,err,k,y=newton(f,df,p0,delta,max1)% f 是非线性系数% df 是 f 的微商% p0 是初始值% dalta 是给定允许误差% max1 是迭代的最大次数% p1 是牛顿法求得的方程近似解% err 是 p0 误差估计% k 是迭代次数p0
4、,feval(f,p0)for k=1:max1- 3 -p1=p0-feval(f,p0)/feval(df,p0);err=abs(p1-p0);p0=p1;p1,err,k,y=feval(f,p1)if(err 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- 7 - 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.50
5、00 -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 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.400
6、0 0.3000 0.2000 0.1000 0追赶法为一种特殊的 LU 分解法。追赶法是求解三对角矩阵的常用方法，但从整体编程角度分析，其程序编写较迭代法复杂，但通用性较好。追赶法求解三对角矩阵不但节省存储单元，而且可以减少计算量，是工程技术中比较常用的数学工具。- 8 -三、数值积分1、题目卫星轨道是一个椭圆，椭圆周长的计算公式是 , dacS202sin)(14这里 是椭圆的半长轴, 是地球中心与轨道中心(椭圆中心)的距离, 记 为近地ac h点距离, 为远地点距离, 公里为地球半径, 则 , H6371R,22RHc某人造卫星近地点距离 公里, 远地点距离 公里, 试用 Romberg
7、5h483方法求卫星轨道的周长，精确到 。602.数学原理龙贝格方法是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。 作为一种外推算法, 它在不增加计算量的前提下提高了误差的精度。龙贝格方法的主要过程是将粗糙的梯形公式 逐步加工成精度较高的辛)(fTn普森公式 和科特斯公式 的方法称为龙贝格方法。)(fSn )(fCn复化梯形公式 )()2)()2)( 110101 iinikknkk xffhxffxfT在复化梯形公式中，每个内节点 既是前一个小区间的终点，又,121n是后一个小区间的起点，因此上式可以改写为 1)()(2)(nkn bfxfafhfT复
8、化梯形公式余项)(12“fhabE,ba- 9 -复化梯形公式的递推公式为 )(21bfabT复化辛普森求积公式 )(1022iiinxfh与复化梯形公式 类似，每个内节点 需用两次，因此有)(fTn n,10211 )(4)(2)(6knkn bfxfxfafhfS显然复化辛普森公式在 n 趋于无穷大时，他的收敛速度比复化梯形公式更快。以 表示二分 k 次后求得的梯形值，且以 ()kmT表示序列 ()0kT的 m 次加速()0kT度，理查森外推法的递推公式可写成 ()(1)()14,2mkkkmT龙贝格算法的计算过程如下：(1)取 求0,khba(0)()2hTfab(2)利用变步长梯形公式
9、 ()0k，其中 k 为区间的二分次数，即 12()(1)00 102()2kkk jjbaTfffx或 12()(1)000(1)22kkk kjbabaTfffj(3)依横行次序求加速值，逐个求出的第 k 行其余各元素 (),)jTk(4)当相邻对角元素之差的绝对值小于预先给定的精度时，终止计算。- 10 -表 3-1 龙贝格算法递推表k h )(0kT)(1k)(2kT)(3k)(4kT0 b-a )(1 2ab)1(0)0(12 4)2(T)()0(2T3 8)3(0)2(1)1()0(34 16ab)4()3()2()1()0(4T3.程序设计function R=romberg(f
10、,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);endR(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);- 11 -endend4.结果分析与讨论本题根据算法原理在 matlab 中编写完龙贝格算法的自定义程序后，直接输入符合格式的函数积分就可得到相应轨道周长。调用MATLAB龙贝格算法的函数后可算得R = romberg(4*7800*sqrt(1-(973.5/7880)2*sin(x)2),0,pi/2,6)计算出来得出R=49136.836545由此可得精10 -6确到的卫星轨道周长约为49136.836545公里。通过本次编程，我对龙贝格算法的公式和步骤有了进一步的掌握，知道了使用龙贝格计算积分是十分方便的，知道 ，就可以知道积分值是多少了，()fx并且误差在误差范围之内，这在数学的计算中是十分重要的，从而解决了许多实际的工程问题。
