數(shù)值計(jì)算大作業(yè)——劉峰

1、課 程 設(shè) 計(jì)課程名稱：數(shù)值分析設(shè)計(jì)題目：數(shù)值計(jì)算大作業(yè)學(xué) 號(hào)：S315070064姓 名：劉峰完成時(shí)間：2015年10月25日題目一、非線性方程求根1.題目 假設(shè)人口隨時(shí)間和當(dāng)時(shí)人口數(shù)目成比例連續(xù)增長，在此假設(shè)下人口在短期內(nèi)的增長建立數(shù)學(xué)模型。（1）如果令表示在時(shí)刻的人口數(shù)目，表示固定的人口出生率，則人口數(shù)目滿足微分方程，此方程的解為；（2）如果允許移民移入且速率為恒定的，則微分方程變成，此方程的解為；假設(shè)某地區(qū)初始有1000000人，在第一年有435000人移入，又假設(shè)在第一年年底該地區(qū)人口數(shù)量1564000人，試通過下面的方程確定人口出生率，精確到；且通過這個(gè)數(shù)值來預(yù)測(cè)第二年年末的人口數(shù)

2、，假設(shè)移民速度保持不變。2.數(shù)學(xué)原理采用牛頓迭代法，牛頓迭代法的數(shù)學(xué)原理是，對(duì)于方程，如果是線性函數(shù)，則它的求根是很容易的，牛頓迭代法實(shí)質(zhì)上是一種線性化方法，其基本思想是將非線性方程逐步歸結(jié)為某種線性方程來求解。設(shè)已知方程有近似根(假定)，將函數(shù)在點(diǎn)進(jìn)行泰勒展開，有于是方程可近似地表示為這是個(gè)線性方程，記其根為，則的計(jì)算公式為，這就是牛頓迭代法，簡稱牛頓法。3.程序設(shè)計(jì)作出函數(shù)的圖像，大概估計(jì)出根的位置fplot('1000*exp(x)+(435*x)*(exp(x)-1)-1564',0 3);grid大概估計(jì)出初始值x=0.5function p1,err,k,y=new

3、ton(f,df,p0,delta,max1)% f是非線性系數(shù)% df是f的微商% p0是初始值% dalta是給定允許誤差% max1是迭代的最大次數(shù)% p1是牛頓法求得的方程近似解% err是p0誤差估計(jì)% k是迭代次數(shù)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.結(jié)果分析與討論在MATLAB中的command window輸入newton('f','df',1.2,10(-4),10)運(yùn)行后得出結(jié)果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運(yùn)算后的結(jié)果為，通過這個(gè)數(shù)值來預(yù)測(cè)第二年年末的人口數(shù)，t=2時(shí)候?qū)τ趯?shí)踐表明，當(dāng)初始值難以確定時(shí)，迭代法就不一定收斂了，因此要根據(jù)問題實(shí)際背景或者二分法先得一個(gè)較好的初始值，然后再進(jìn)行迭代；再者迭代函數(shù)選擇不合適的話，采用不動(dòng)點(diǎn)迭代法也有可能出現(xiàn)不收斂的情況；因此我采用的是牛頓法。題目二：線性方程組求解1.題目 假設(shè)一個(gè)物體可以位于個(gè)等距

6、點(diǎn)的任意位置，當(dāng)物體在位置時(shí)，它只能等可能的移動(dòng)到或者，而不能直接移動(dòng)到其他任何位置，概率表示物體從位置開始在到達(dá)右端點(diǎn)之前到達(dá)左端點(diǎn)的概率，顯然，且有既有下面方程組：取對(duì)方程組進(jìn)行求解（迭代法或者直接法）。2.數(shù)學(xué)原理 在解微分方程的邊值問題、熱傳導(dǎo)方程以及船體數(shù)學(xué)放樣中建立的三次樣條函數(shù)等工程技術(shù)問題時(shí)，經(jīng)常遇到下面形式的線性方程組： =方程簡記，該線性方程稱為三對(duì)角線方程組，其系數(shù)矩陣A滿足條件所以為弱對(duì)角陣可以采用追趕法進(jìn)行計(jì)算，利用三對(duì)角矩陣的LU分解建立計(jì)算量更少的線性方程組求解公式。將系數(shù)矩陣A進(jìn)行克勞特分解，即A分解為下三角矩陣和單位上三角矩陣的乘積；A= 其中，為待定系數(shù)，直

7、接利用矩陣乘法公式可得，于是推得計(jì)算，的公式，；，；，；由此計(jì)算出L和U中的全部元素，完成了系數(shù)矩陣A的克勞特分解。求解線性方程組等價(jià)于求解和。因而得到解三對(duì)角線性方程組的追趕法公式(1) 計(jì)算的遞推公式：(2) 解(3) 解我們將計(jì)算系數(shù)和稱為追的過程，將計(jì)算方程組的解稱為趕的過程。整個(gè)過程為追趕法的思想。3.程序設(shè)計(jì)function x=chase (a,b,c,f) %求解線性方程組Ax=f，其中A是三對(duì)角陣 %a 是矩陣A的下對(duì)角線元素 a（1）=0 %b 是矩陣A的對(duì)角線元素 %c 是矩陣A的上對(duì)角線元素 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.結(jié)果分析與討論A的系數(shù)矩陣為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中調(diào)用之前保存的迭代法函數(shù)function，在命令窗口中輸入：chase(a,b,c,f)回車得到結(jié)果：>> 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分解法。追趕法是求解三對(duì)角矩陣的常用方法，但從整體編程角度分析，其程序編寫較迭代法復(fù)雜，但通用性較好。追趕法求解三對(duì)角矩陣不但節(jié)省存儲(chǔ)單元，而且可以減少計(jì)算量，是工程技術(shù)中比較常用的數(shù)學(xué)工具。三、數(shù)值積分1、題目

12、衛(wèi)星軌道是一個(gè)橢圓，橢圓周長的計(jì)算公式是, 這里是橢圓的半長軸, 是地球中心與軌道中心(橢圓中心)的距離, 記為近地點(diǎn)距離, 為遠(yuǎn)地點(diǎn)距離, 公里為地球半徑,則, 某人造衛(wèi)星近地點(diǎn)距離公里,遠(yuǎn)地點(diǎn)距離公里, 試用Romberg方法求衛(wèi)星軌道的周長，精確到。2.數(shù)學(xué)原理龍貝格方法是在梯形公式、辛普森公式和柯特斯公式之間的關(guān)系的基礎(chǔ)上，構(gòu)造出一種加速計(jì)算積分的方法。 作為一種外推算法, 它在不增加計(jì)算量的前提下提高了誤差的精度。 龍貝格方法的主要過程是將粗糙的梯形公式逐步加工成精度較高的辛普森公式和科特斯公式的方法稱為龍貝格方法。復(fù)化梯形公式 在復(fù)化梯形公式中，每個(gè)內(nèi)節(jié)點(diǎn)既是前一個(gè)小區(qū)間的終點(diǎn)，又

13、是后一個(gè)小區(qū)間的起點(diǎn)，因此上式可以改寫為復(fù)化梯形公式余項(xiàng) 復(fù)化梯形公式的遞推公式為復(fù)化辛普森求積公式與復(fù)化梯形公式類似，每個(gè)內(nèi)節(jié)點(diǎn)需用兩次，因此有顯然復(fù)化辛普森公式在n趨于無窮大時(shí)，他的收斂速度比復(fù)化梯形公式更快。以表示二分k次后求得的梯形值，且以表示序列的m次加速度，理查森外推法的遞推公式可寫成龍貝格算法的計(jì)算過程如下：(1) 取求(2) 利用變步長梯形公式，其中k為區(qū)間的二分次數(shù)，即或(3) 依橫行次序求加速值，逐個(gè)求出的第k行其余各元素(4) 當(dāng)相鄰對(duì)角元素之差的絕對(duì)值小于預(yù)先給定的精度時(shí)，終止計(jì)算。表3-1龍貝格算法遞推表kh0b-a12343.程序設(shè)計(jì)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.結(jié)果分析與討論本題根據(jù)算法原理在matlab中編寫完龍貝格算法的自定義程序后，直接輸入符合格式的函數(shù)積分就可得到相應(yīng)軌道周長。調(diào)用MATLAB龍貝格算法的函數(shù)后可算得 R = romber