数值计算大作业刘

课程设计课程名称：数值分析设计题目：数值计算大作业学号：S315070064姓名：刘峰完成时间：2015年10月25日-1-题目一、非线性方程求根1.题目假设人口随时间和当时人口数目成比例连续增长，在此假设下人口在短期内的增长建立数学模型。（1）如果令()Nt表示在t时刻的人口数目，表示固定的人口出生率，则人口数目满足微分方程()()dNtNtdt，此方程的解为0()=tNtNe；（2）如果允许移民移入且速率为恒定的v，则微分方程变成()()dNtNtvdt，此方程的解为0()=+(1)ttvNtNee；假设某地区初始有1000000人，在第一年有435000人移入，又假设在第一年年底该地区人口数量1564000人，试通过下面的方程确定人口出生率，精确到410；且通过这个数值来预测第二年年末的人口数，假设移民速度v保持不变。4350001564000=1000000(1)ee2.数学原理采用牛顿迭代法，牛顿迭代法的数学原理是，对于方程0)(xf，如果)(xf是线性函数，则它的求根是很容易的，牛顿迭代法实质上是一种线性化方法，其基本思想是将非线性方程0)(xf逐步归结为某种线性方程来求解。设已知方程0)(xf有近似根kx(假定0)(xf)，将函数)(xf在点kx进行泰勒展开，有.))(()()(kkkxxxfxfxf于是方程0)(xf可近似地表示为0))(()(kkxxxfxf这是个线性方程，记其根为1kx，则1kx的计算公式为)()(1kkkkxfxfxx，,,2,1,0k-2-这就是牛顿迭代法，简称牛顿法。3.程序设计作出函数的图像，大概估计出根的位置fplot('1000*exp(x)+(435*x)*(exp(x)-1)-1564',[03]);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,feval('f',p0)fork=1:max1p1=p0-feval('f',p0)/feval('df',p0);-3-err=abs(p1-p0);p0=p1;p1,err,k,y=feval('f',p1)if(errdelta)|(y==0),break,endp1,err,k,y=feval('f',p1)endfunctiony=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('f','df',1.2,10^(-4),10)运行后得出结果p0=0.5000p1=0.1679err=0.3321k=1y=9.2415e+004p1=0.1031err=0.0648k=2y=2.7701e+003p1=0.1010err=0.0021k=3y=2.6953p1=0.1010err=2.0129e-006k=4y=2.5576e-006ans=0.1010运算后的结果为1010.0，通过这个数值来预测第二年年末的人口数，0.10100.1010435000f(t)=1000000(1)0.1010tteet=2时候对于f()2187945.865x实践表明，当初始值难以确定时，迭代法就不一定收敛了，因此要根据问题实际背景或者二分法先得一个较好的初始值，然后再进行迭代；再者迭代函数选择不合适的话，采用不动点迭代法也有可能出现不收敛的情况；因此我采用的是牛顿法。-4-题目二：线性方程组求解1.题目假设一个物体可以位于1n个等距点01,,,nxxx的任意位置，当物体在ix位置时，它只能等可能的移动到1ix或者+1ix，而不能直接移动到其他任何位置，概率ip表示物体从位置ix开始在到达右端点nx之前到达左端点0x的概率，显然01,0npp，且有-1+111=+1,2,,122iiipppin，既有下面方程组：1211100211112221110221101221012nppp取10n对方程组进行求解（迭代法或者直接法）。2.数学原理在解微分方程的边值问题、热传导方程以及船体数学放样中建立的三次样条函数等工程技术问题时，经常遇到下面形式的线性方程组：nnnnnbacbabacb1112211nxxxxn121=nndddd121方程简记Axd，该线性方程称为三对角线方程组，其系数矩阵A满足条件110,,,0,2,,10iiiiinnbcbacacinbc-5-所以为弱对角阵可以采用追赶法进行计算，利用三对角矩阵的LU分解建立计算量更少的线性方程组求解公式。将系数矩阵A进行克劳特分解，即A分解为下三角矩阵和单位上三角矩阵的乘积；A=nnnnnbacbabacb1112211=nnnn112211111121n其中i，i，i为待定系数，直接利用矩阵乘法公式可得11b，111c，iia，iiiib1，,,,3,2niiiic，,1,,3,2ni于是推得计算i，i，i的公式11b，111/bc；ii，1iiiib，,,,3,2ni；iiic/，ni,,3,2；由此计算出L和U中的全部元素，完成了系数矩阵A的克劳特分解。求解线性方程组dAx等价于求解dLy和yUx。因而得到解三对角线性方程组的追赶法公式(1)计算ib的递推公式：1111,,2,3,1iiiiicbcbinbbab(2)解Lyd11111,,2,3,iiiiiiiydbydaybinab(3)解Uxy1,,1,,1nniiiixyxyxinb我们将计算系数121n和121nyyy称为追的过程，将计算方程组的解11nnxxx称为赶的过程。整个过程为追赶法的思想。-6-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:1a(i+1)=a(i);endendc(1)=c(1)/b(1);f(1)=f(1)/b(1);fori=2:n-1b(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));fori=n-1:-1:1f(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,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]-7-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.5000b=1111111111c=-0.5000-0.5000-0.5000-0.5000-0.5000-0.5000-0.5000-0.5000-0.50000f=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分解法。追赶法是求解三对角矩阵的常用方法，但从整体编程角度分析，其程序编写较迭代法复杂，但通用性较好。追赶法求解三对角矩阵不但节省存储单元，而且可以减少计算量，是工程技术中比较常用的数学工具。-8-三、数值积分1、题目卫星轨道是一个椭圆，椭圆周长的计算公式是dacaS2022sin)(14,这里a是椭圆的半长轴,c是地球中心与轨道中心(椭圆中心)的距离,记h为近地点距离,H为远地点距离,6371R公里为地球半径,则2,22RHhHhac,某人造卫星近地点距离536h公里,远地点距离2483H公里,试用Romberg方法求卫星轨道的周长，精确到610。2.数学原理龙贝格方法是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。作为一种外推算法,它在不增加计算量的前提下提高了误差的精度。龙贝格方法的主要过程是将粗糙的梯形公式)(fTn逐步加工成精度较高的辛普森公式)(fSn和科特斯公式)(fCn的方法称为龙贝格方法。复化梯形公式)]()([2)()(2)(1101101iinikknkkknxfxfhxfxfxxfT在复化梯形公式中，每个内节点,,,121nxxx既是前一个小区间的终点，又是后一个小区间的起点，因此上式可以改写为11)()(2)(2)(nkknbfxfafhfT复化梯形公式余项)(122fhabE],[ba复化梯形公式的递推公式为)]()([21bfafabT-9-复化辛普森求积公式)(22110212iiinnxfhTT与复化梯形公式)(fTn类似，每个内节点nxxx,,,21需用两次，因此有102111)()(4)(2)(6)(nkknkknbfxfxfafhfS显然复化辛普森公式在n趋于无穷大时，他的收敛速度比复化梯形公式更快。以()0kT表示二分k次后求得的梯形值，且以()kmT表示序列()0kT的m次加速度，理查森外推法的递推公式可写成()(1)()1141,1,2,4141mkkkmmmmmTTTk龙贝格算法的计算过程如下：(1)取0,,khba求(0)0()()2hTfafb(2)利用变步长梯