数值计算大作业

数值计算大作业题目一、非线性方程求根1.题目假设人口随时间和当时人口数目成比例连续增长，在此假设下人口在短期内的增长建立数学模型。

（1）如果令()N t 表示在t 时刻的人口数目，β表示固定的人口出生率，则人口数目满足微分方程()()dN t N t dt β=，此方程的解为0()=tN t N e β； （2）如果允许移民移入且速率为恒定的v ，则微分方程变成()()dN t N t vdt β=+， 此方程的解为0()=+(1)t t vN t N e e βββ-；假设某地区初始有1000000人，在第一年有435000人移入，又假设在第一年年底该地区人口数量1564000人，试通过下面的方程确定人口出生率β，精确到410-；且通过这个数值来预测第二年年末的人口数，假设移民速度v 保持不变。

4350001564000=1000000(1)e e βββ+-2.数学原理采用牛顿迭代法，牛顿迭代法的数学原理是，对于方程0)(=x f ，如果)(x f 是线性函数，则它的求根是很容易的，牛顿迭代法实质上是一种线性化方法，其基本思想是将非线性方程0)(=x f 逐步归结为某种线性方程来求解。

设已知方程0)(=x f 有近似根k x (假定0)(≠'x f )，将函数)(x f 在点k x进行泰勒展开，有.))(()()(⋅⋅⋅+-'+≈k k k x x x f x f x f于是方程0)(=x f 可近似地表示为))(()(=-'+k k x x x f x f这是个线性方程，记其根为1k x +，则1k x +的计算公式为)()(1k k k k x f x f x x '-==+，,,2,1,0⋅⋅⋅=k这就是牛顿迭代法，简称牛顿法。

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,feval('f',p0) for k=1:max1p1=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,endp1,err,k,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)/x^2);4.结果分析与讨论newton('f','df',1.2,10^(-4),10) 运行后得出结果 p0 =0.5000p1 =0.1679 err =0.3321 k =1 y =9.2415e+004 p1 =0.1031 err =0.0648 k =2 y =2.7701e+003 p1 =0.1010 err =0.0021 k =3 y =2.6953p1 =0.1010 err =2.0129e-006 k =4 y = 2.5576e-006 ans =0.1010运算后的结果为1010.0=β，通过这个数值来预测第二年年末的人口数，0.10100.1010435000f(t)=1000000(1)0.1010t te e +-t=2时候对于f ()2187945.865x =实践表明，当初始值难以确定时，迭代法就不一定收敛了，因此要根据问题实际背景或者二分法先得一个较好的初始值，然后再进行迭代；再者迭代函数选择不合适的话，采用不动点迭代法也有可能出现不收敛的情况；因此我采用的是牛顿法。

题目二：线性方程组求解1.题目假设一个物体可以位于1n +个等距点01,,,n x x x 的任意位置，当物体在i x位置时，它只能等可能的移动到1i x -或者+1i x ，而不能直接移动到其他任何位置，概率i p 表示物体从位置ix开始在到达右端点n x 之前到达左端点0x 的概率，显然01,0n p p ==，且有-1+111=+1,2,,122i i i p p p i n =-， 既有下面方程组：1211100211112221110221101221012n p p p-⎡⎤-⎢⎥⎢⎥⎢⎥--⎢⎥⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦--⎢⎥⎢⎥⎢⎥-⎢⎥⎣⎦取10n =对方程组进行求解（迭代法或者直接法）。

2.数学原理在解微分方程的边值问题、热传导方程以及船体数学放样中建立的三次样条函数等工程技术问题时，经常遇到下面形式的线性方程组：⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡---n n n n n b a c b a b a c b 1112211 ⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-n x x x x n 121 =⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-n n d d d d 121 方程简记A x d =，该线性方程称为三对角线方程组，其系数矩阵A 满足条件110,,,0,2,,1i i i i i n n b c b a c a c i n b c >>≥+≠=->>所以为弱对角阵可以采用追赶法进行计算，利用三对角矩阵的LU 分解建立计算量更少的线性方程组求解公式。

将系数矩阵A 进行克劳特分解，即A 分解为下三角矩阵和单位上三角矩阵的乘积；A=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡---n nn n n b a c b a b a c b 1112211=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡--n n n n αγαγαγα11221 ⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-1111121n βββ其中i α，i β，i γ为待定系数，直接利用矩阵乘法公式可得11α=b ，111βα=c ， ii a γ=，ii i i b αβγ+=-1，,,,3,2n i ⋅⋅⋅=ii i c βα=，,1,,3,2-⋅⋅⋅=n i于是推得计算i α，i β，i γ的公式11b =α，111/b c =β；i i αγ=，1--=i i i i b βαα，,,,3,2n i =；i i i c αβ/=，n i ,,3,2 =；由此计算出L 和U 中的全部元素，完成了系数矩阵A 的克劳特分解。

求解线性方程组d Ax =等价于求解d Ly =和y Ux =。

因而得到解三对角线性方程组的追赶法公式 (1)计算i 的递推公式：()1111,,2,3,1i i i i i c b c b i n -==-=-(2)解Ly d =())11111,,2,3,i i i i i i i y d b y d a y b i n --==--=(3)解Ux y = 1,,1,,1n n i i ii x y x y x i n +==-=-我们将计算系数121n βββ-→→→和121n y y y -→→→称为追的过程，将计算方程组的解11n n x x x -→→→称为赶的过程。

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

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==length(a)for i=n-1:-1:1a(i+1)=a(i);endendc(1)=c(1)/b(1);f(1)=f(1)/b(1);for i=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));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,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.5000b =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.4000 0.3000 0.20000.1000 0追赶法为一种特殊的LU分解法。