《數(shù)值分析》上機實驗報告

1、數(shù)值分析上機實驗報告?數(shù)值分析?上機實驗報告1. 用Newton法求方程X7-X4+14=0在(0.1, 1.9)中的近似根(初始近似值取為區(qū)間端點，迭代6次或誤差小于0.00001)。1.1理論依據(jù)：設(shè)函數(shù)在有限區(qū)間a, b上二階導(dǎo)數(shù)存在，且滿足條件1. f(x)f(b) 02. f(x)在區(qū)間a, b上不變號3f(x)0;4f (c) | f .(x) |,其中 c是a,b中使 mir(| f .(a), f .(b) |)到達(dá)的一個 b a那么對任意初始近似值x0a,b,由Newton迭代過程f (xk)Xk1(Xk) Xk k ,k 0,1,2,3f'(Xk)所生的迭代序列Xk

2、平方收斂于方程f(x) 0在區(qū)間a,b上的惟一解令74f (x) x 28x14, f (0.1)0, f (1.9)0f (x) 7x6112x3 7x3(x316)0f (x)42x5336x242x2(x38)0f (1.9) f (1.9)0故以1.9為起點f(Xk)Xk 1 Xkf (Xk)x0 1.9如此一次一次的迭代，逼近X的真實根。當(dāng)前后兩個的差=£時， 就認(rèn)為求出了近似的根。本程序用 Newt on法求代數(shù)方程(最高次數(shù)不 大于10)在(a,b)區(qū)間的根。1.2 C語言程序原代碼:#i nclude<stdio.h>#in clude<math.h

3、> mai n()double x2,f,f1;double x1=1.9;/ 取初值為 1.9dox2=x1;f=pow(x2,7)-28*pow(x2,4)+14; f1=7*pow(x2,6)-4*28*pow(x2,3);/限制循環(huán)次數(shù)x仁 x2-f/f1;while(fabs(x1-x2)>=0.00001|x1<0.1); printf("計算結(jié)果:x=%fn",x1);1.3運行結(jié)果：1.4 MATLAB上機程序fun cti on y=Newt on( f,df,x0,eps,M) d=0;for k=1:Mif feval(df,x0)=

4、0 d=2;breakelse x1=x0-feval(f,x0)/feval(df,x0);ende=abs(x1-x0);x0=x1;if e<=eps&&abs(feval(f,x1)v=eps d=1;breakend endif d=1y=x1;elseif d=0y='迭代M次失敗'；elsey= '奇異 'endfunction y=df(x)y=7*xA6-28*4*xA3;Endfunction y=f(x)y=xA7-28*xA4+14;End>> x0=1.9;>> eps=0.00001;&g

5、t;> M=100;>> x=Newton('f','df',x0,eps,M);>> vpa(x,7)1.5 問題討論：1. 使用此方法求方解，用誤差來控制循環(huán)迭代次數(shù)，可以在誤差允許的范圍 內(nèi)得到比較理想的計算結(jié)果。 此程序的缺乏之處是， 所要求解的方程必須滿足上 述定理的四個條件，但是第二和第四個條件在計算機上比較難以實現(xiàn)。2. Newton迭代法是一個二階收斂迭代式，他的幾何意義Xi+1是Xi的切線與 x軸的交點，故也稱為切線法。它是平方收斂的，但它是局部收斂的，即要求初始 值與方程的根充分接近，所以在計算過程中需要先確定初

6、始值。3. 此題在理論依據(jù)局部，討論了區(qū)間(0.1,1.9)兩端點是否能作為Newton迭代 的初值，結(jié)果發(fā)現(xiàn) 0.1 不滿足條件，而 1 . 9滿足，能作為初值。另外，該程序簡 單，只有一個循環(huán)，且為順序結(jié)構(gòu)，故采用 do-while 循環(huán)。當(dāng)然也可以選擇 for 和 while 循環(huán)。2. 函數(shù)值如下表:X12345f(x)o0.693147181.o9861231.38629441.6094378X67891of(x)1.79175951.94591012.0794452.19722462.3025851f(x)f (1)=1f (10)=0.1試用三次樣條插值求f(4.563)及f (

7、4.563)的近似值2.1理論依據(jù)S(x) Mj(Xj x)31 _6hj 1(x Xj 1)3Mj L- (yj1 Mj6hj 1.2h, 1 X, x1士)(+ )6 hj 1(力.2hj 1 x Xj 1 叫士)(十-)6 hj 1這里hj 1XjXj，所以只要求出就能得出插值函數(shù)S(x)。MoM1doa這里的方法為:dNdodj_6( _hoho6 (yj1 yj hj 1 hjhj61h yN h (yN h N 1hN 1hj 1hj 1 hjyoyo)yjhj 1yN 1 )丄)(j 1,2,L ,Ni)最終歸結(jié)為求解一個三對角陣的解。用追趕法解三對角陣的方法如下:idi, li

8、 ii 1bi 1li iC ,AbiG82b2OC2O3n 1O1l11l2 1O1Oln 112 2O On 1n 1nLU0 1 canbnn 1Ld卄1LUxd,即,假設(shè)記M，那么由Ld得Uxn11d111X11l21MMOOMMO OMMn 1n 1MMln1ndnnXnn綜上可得求解方程Ax=d的算法:i 1 di 1 1i 1 i i 1,2,3,L ,n 1Xn, Xi丄沁，i n 1,L ,2,1ni2.2 C語言程序代碼:#i nclude<stdio.h>#in clude<math.h> void mai n()int i,j,m, n,k,p;

9、double q10,p10,s4,g4,x0,x1,g0=1,g9=0.1;double s1010;double a10,b10,c10,d10,e10,x10,h9,u9,r9;double f10=0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.1972246,2.3025851; printf("請依次輸入 xi:n");for(i=0;i<=9;i+)scanf("%lf",&ei);/求 h 矩陣for(n=0;n<=8;

10、n+)h n=e n+1-e n;d0=6*(f1-f0)/h0-g0)/h0;d9=6*(g9-(f9-f8)/h8)/h8;for(j=0;j<=7;j+) dj+1=6*(fj+2-fj+1)/hj+1-(fj+1-fj)/hj)/(hj+hj+1);for(m=1;m<=8;m+) um=hm-1/(hm-1+hm);for(k=1;k<=8;k+)rk=hk/(hk-1+hk);for(i=0;i<=9;i+) / 求 u 矩陣 for(p=0;p<=9;p+)sip=0;if(i=p)sip=2;s01=1;s98=1;for(i=1;i<=8;

11、i+)sii-1=ui;sii+1=ri;printf(" 三對角矩陣為 :n");for(i=0;i<=9;i+)for(p=0;p<=9;p+)/求 r 矩陣 printf("%5.2lf",sip);if(p=9)printf("n");printf(" 根據(jù)追趕法解三對角矩陣得： n");a0=s00;b0=d0;for(i=1;i<9;i+)ci=sii-1/ai-1;/求 d 矩陣ai=sii-si-1i*ci;bi=di-ci*bi-1;if(i=8)p10=bi;q10=ai; x

12、9=p10/q10;printf("M10=%lfn",x9);for(i=9;i>=1;i-) xi-1=(bi-1-si-1i*xi)/ai-1;printf("M%d=%lfn",i,xi-1);printf("可得s(x)在區(qū)間4,5上的表達(dá)式；n");printf("將 x=4.563 代入得：n");x0=5-4.563;x1=4.563-4;s4=x3*pow(x0,3)/6+x4*pow(x1,3)/6+(f3-x3/6)*(5-4.563)+(f4-x4/6)*(4.563 -4);g4=-

13、x3*pow(x0,2)/2+x4*pow(x1,2)/2-(f3-x3/6)+(f4-x4/6);printf("計算結(jié)果：f(4.563)的函數(shù)值是：lfnf(4.563)的導(dǎo)數(shù)值是：lfn",s4,g4);2.3 運行結(jié)果：2.4 問題討論1. 三次樣條插值效果比Lagrange插值好，沒有Runge現(xiàn)象，光滑性較好。2. 此題的對任意劃分的三彎矩插值法可以解決非等距節(jié)點的一般性問題。3. 編程過程中由于定義的數(shù)組比較多，需要仔細(xì)弄清楚各數(shù)組所代表的參 數(shù)，要注意各下標(biāo)代表的含義，特別是在用追趕法計算的過程中 。33.用 Romberg算法求 3xx1'4(5

14、x 7)sinx2dx(允許誤差0.00001).3.1理論依據(jù)：Romberg算法的計算步驟如下:(1) 先求出按梯形公式所得的積分值f(a)f(b)(2) 把區(qū)間2等分，求出兩個小梯形面積之和，記為 T,，即久晉f(a) f(b) 2f (穿)這樣由外推法可得T2(0) , t2(0)(1)扌聲1 (；)24T, T,(0)n(3) 把區(qū)間再等分(即22等分)，得復(fù)化梯形公式T,，由T,與T外推可(2) (1) 2 (1) (0)得T2 紅 匚，T3(0) 4 T22 T2 ，如此，假設(shè)已算出2k等分的復(fù)化梯形公式4 141TV，那么由Richards on外推法，構(gòu)造新序列4mT (k)

15、 Tk 1)Tnmk11)4 Tmm Tm ，m=1,2,1, k=1,2,l-m+141最后求得Tl(0)(4) T|(0)丁黑或 |T|(0)Ti(01 |就停止計算，否那么回到(3)，計算T1(l 1)，般可用如下算法:T1(0)T1(l)(k 1)T m 1b a f(a)21 (i 1) b 尹 P i14mT,k) T,k1)f(b)211afa(2i1,2,L ,l , k 1,2,L ,l m 1其具體流程如下，并全部存入第一列通常計算時，用固定l=N來計算，一般1=4或5即能到達(dá)要求。3.2 C語言程序代碼：#in clude<math.h>#i nclude&l

16、t;stdio.h>double f(double x)/計算 f(x)的值double z;z=pow(3,x)*pow(x,1.4)*(5*x+7)*si n(x*x);return(z);mai n() double t2020,s,e=0.00001,a=1,b=3;下為的差是int i,j,l,k;t01=(b-a)*(f(b)+f(a)/2;/romberg 算法t11=(b-a)*(f(b)+2*f(b+a)/2)+f(a)/4; t02=(a*t11-t01)/(4-1);j=3;for(l=2;fabs(t0j-1-t0j-2)>=e;l+)for(k=1,s=0

17、;k<=pow(2,l-1);k+)s+=f(a+(2*k-1)*(b-a)/pow(2,l);/判斷前后兩次所得的T(0)否符合要求，如果符合精度要求那么停止循環(huán)tl1=(tl-11+(b-a)*s/pow(2,l-1)/2;for(i=l-1,j=2;i>=0;i-,j+) tij=(pow(4,j-1)*ti+1j-1-tij-1)/(pow(4,j-1)-1);if(t01<e)prin tf("t=%0.6fn",t01);elseprintf(" 用 Romberg算法計 算函數(shù) 所得近 似結(jié)果 nf(x)=%0.6fn",

18、t0j-1);3.3運行結(jié)果:'D WC+ -EXERCISE3XDtbug3.eRThe yeisLiil't is "440,5301 Piess bfi¥to continue3.4 MATLAB上機程序function T,n=mromb(f,a,b,eps)if narginv4,eps=1e-6; endh=b-a;R(1,1)=(h/2)*(feval(f,a)+feval(f,b);n=1;J=O;err=1;while (err>eps)J=J+1;h=h/2;S=0;for i=1: nx=a+h*(2*i-1);S=S+feval

19、(f,x);endR(J+1,1)=R(J,1)/2+h*S;for k=1:JR(J+1,k+1)=(4Ak*R(J+1,k)-R(J,k)/(4Ak-1);enderr=abs(R(J+1,J+1)-R(J+1,J);n=2*n;endR；T=R(J+1,J+1);format lo ngf=(x)(3.Ax)*(x.A1.4)*(5*x+7)*si n(x*x);T, n=mromb(f,1,3,1.e-5)3.5問題討論：1. Romberge算法的優(yōu)點是：把積分化為代數(shù)運算，而實際上只需求Ti(i),以后用遞推可得.算法簡單且收斂速度快，一般4或5次即能到達(dá)要求。2. Romberg

20、e 算法的缺點是 :對函數(shù)的光滑性要求較高， 計算新分點 的值時，這些數(shù)值的個數(shù)成倍增加。3. 該程序較為復(fù)雜，涉及函數(shù)定義，有循環(huán)，而且循環(huán)中又有判 斷，編寫時需要注意該判斷條件是處于循環(huán)中， 當(dāng)?shù)竭_(dá)要求時跳出循 環(huán)，終止運算。4. 函數(shù)的定義可放在主函數(shù)前也可在主程序后面。本程序采用的 后置方式。4.用定步長四階Runge-Kutta求解dy1 /dt 1dy2 /dt y3dy3 /dt 1000 1000y2 100y3y1(0) 0y2(0) 0y3(0) 0h=0.0005,打印 yi(0.025) , yi(0.045) , yi(0.085) , yi(0.1),(i=1， 2

21、， 3)4.1 理論依據(jù)：Ru nge_Kutta采用高階單步法，這里不是先按Taylor公式展開，而是先寫成tn處附近的值的線性組合(有待定常數(shù))再按 Taylor公式展開，然后確定待定常數(shù)，這就是 Runge-Kutta 法的思想方法。此題采用四階古典的Run ge-Kutta公式：Yn 1 Yn K1 3K2 3K3 K4/8K1hF ( xn ,Yn )K2hF(xnh/3,YnhK1/3)K3hF(xn2h/3,YnhK1 /3hK2)K4hF(xnh,Yn hK1 hK2hK3)4.2 C 語言程序代碼：#include<stdio.h>void fun(double

22、x4,double y4,double h) y1=1*h;y2=x3*h;y3=(1000-1000*x2-100*x2-100*x3)*h;/微分方程向量函數(shù) void main() double Y54,K54,m,z4,e=0.0005;double y5=0,0.025,0.045,0.085,0.1;int i,j,k;for(i=1;i<=3;i+)Y1i=0;for(i=1;i<=4;i+)for(j=1;jv=3;j+)Kij=O;for(k=1;k<=5;k+)for(m=yk-1;m<=yk;m=m+e)for(i=1;i<=3;i+)zi=

23、Yki;fun (z,K1,e);for(i=1;i<=3;i+)zi=Yki+e*K2i/2;依此求 K1,K2K3 的值fun (z,K2,e); for(i=1;i<=3;i+)zi=Yki+e*K2i/2;fun( z,K3,e);for(i=1;i<=3;i+)zi=Yki+e*K3i;fun (z,K4,e);/求for(i=1;i<=3;i+)Yki=Yki+(K1i+2*K2i+2*K3i+K4i)/6;YiN+1的值if(k!=5)for(i=1;i<=3;i+)Yk+1i=Yki;printf("計算結(jié)果:n");for(i

24、=1;i<5;i+)for(j=1;j<=3;j+)pri ntf("y%d%4.3f=%-10.8f,",j,yi,Yij);if(j=3)prin tf("n");prin tf("n");4.3運行結(jié)果：025000 yl<0.a45>=0.045000 pl<0-085>=0-a8590S yl<S.10G>=S丄司日y2<0.025>-6.149501 y2<0.045>=0.310966 Lf2<B-B85> =0-559342y3<

25、0.025>-S_33575& y3<0.045>=7-536549 p3<a-685=4-946488Ppeis Artii k電y to continue4.4 問題討論：1.定步長四階 Runge-kutta 方法是一種高階單步法法穩(wěn)定，精度較 高，誤差小且程序相對簡單，存儲量少。不必求出起始點的函數(shù)值， 可根據(jù)精度的要求修改步長，不會由于起始點的誤差造成病態(tài)。2.本程序可以通過修改主程序所調(diào)用的函數(shù)中的表達(dá)式來實現(xiàn) 對其它函數(shù)的任意初值條件求微分計算。3. 程序中運用了大量的 for 循環(huán)語句，因為該公式中涉及大量的求 和，且有不同的函數(shù)和對不同的數(shù)值求

26、值，編程稍顯繁瑣。所以編寫 過程中一定要注意各循環(huán)的次數(shù)，以免出錯。5.12.384122.115237-1.0610741.112336- 0.1135840.7187191.7423823.0678132.0317432.11523719.141823- 3.125432- 1.0123452.1897361.563849- 0.7841651.1123483.123124-1.061074- 3.12543215.5679143.1238482.0314541.836742- 1.0567810.336993- 1.0101031.112336- 1.0123453.12384827.1

27、084374.101011- 3.7418562.101023- 0.71828- 0.037585A- 0.1135842.1897362.0314544.10101119.8979180.431637- 3.1112232.1213141.7843170.7187191.5638491.836742- 3.7418560.4316379.789365- 0.103458- 1.1034560.2384171.742382- 0.784165-1.0567812.101023- 3.111223- 0.10345814.7138465 3.123789- 2.2134743.0678131.

28、1123480.336993- 0.718282.121314-1.1034563.12378930.7193344.446782- 2.0317433.123124-1.010103- 0.0375851.7843170.238417- 2.2134744.44678240.00001b(2.187436933.992318- 25.1734170.846716951.784317- 86.6123431.11012304.719345 - 5.6784392 )用列主元消去法求解 Ax=b。5.1 理論依據(jù)：列主元素消元法是在應(yīng)用Gauss消元法的根底上，憑借長期經(jīng)驗 積累提出的， 是線性

29、方程組一般解法， 目的是為防止在消元計算中使 誤差的擴(kuò)大， 甚至嚴(yán)重?fù)p失了有效數(shù)字使數(shù)據(jù)失真， 而在每次初等變 換前對矩陣作恰當(dāng)?shù)恼{(diào)整，以提高Gauss消元法的數(shù)字穩(wěn)定性，進(jìn)而 提高計算所得數(shù)據(jù)的精確度。 即在每主列中取絕對值最大的元素作主 元，再做對應(yīng)的行交換然后消元求解的方法。具體做法如下：將方陣A和向量b寫成C=( A，b)。將C的第1列中第1行的 元素與其下面的此列的元素逐一進(jìn)行比較，找到最大的元素Cji，將第 j 行的元素與第 1 行的元素進(jìn)行交換，然后通過行變換，將第 1 列中 第2到第n個元素都消成0。將變換后的矩陣c的第二列中第二行 的元素與其下面的此列的元素逐一進(jìn)行比較，找到

30、最大的元素ck(12) ，將第k行的元素與第2行的元素進(jìn)行交換，然后通過行變換，將第 2 列中第3到第n個元素都消成0。以此方法將矩陣的左下局部全都消成 0 后再求解。最終形式如下：a(n) a11 0LOLa(n)a1nMg1M(A,b )MOOMM0L0a(n) nngn5.2 C 語言程序代碼(1)比較該列的元素的絕對值的大小，將絕對值最大的元素通過行 變換使其位于主對角線上；(2)進(jìn)行高斯消去法變換，把系數(shù)矩陣化成上三角形，然后回代求#include "math.h"#include "stdio.h"void Householder(doubl

31、e A99);void expunction(double A99,double b9,double x9);void main()double A99=12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743,2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124,-1.061074,-3.125432,15.567914,3.123848,2.031454,1.8367

32、42,-1.056781,0.336993,-1.010103,1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0. 037585,-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317,0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417,1.742382,-0.78

33、4165,-1.056781,2.101023,-3.111223,-0.103458,14.713847,3.123789,- 2.213474,3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.4 46782,-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,4 0.00001;double b9= 2.1874369,33.992318,-25.173417,0.84671695,1.7

34、84317,-86.612343,1.1101230,4.71 9345,-5.6784392;double x9=0.0;int i,j;Householder(A);printf("n The Results of X are:n");expunction(A,b,x);for(i=1;i<10;i+)printf("X%1d=%fn",i,xi-1);void Householder(double A99)double q9,u9,y9,s,a,kr;int i,j,k;for(i=0;i<7;i+)s=0;for(j=i+1;j<

35、;9;j+) s+=Aji*Aji;s=sqrt(s);a=s*s+fabs(Ai+1i)*s;for(j=0;j<9;j+)if(j<=i) uj=0;else if(j=i+1) uj=Aji+Aji/fabs(Aji)*s;else if(j>i+1) uj=Aji;for(k=0;k<9;k+)yk=0;for(j=0;j<9;j+)yk+=Akj*uj;yk/=a;kr=0;for(k=0;k<9;k+)kr+=yk*uk;kr/=2*a;for(k=0;k<9;k+)qk=yk-kr*uk;for(k=0;k<9;k+)for(j=0

36、;j<9;j+)Akj-=uk*qj+uj*qk;void expunction(double A99,double b9,double x9)int i,j,k;double B910;double z3;double t1=0,t2=0,t3=0;for(i=0;i<8;i+)if(Ai+1i>Aii)for(j=i,k=0;j<i+3;j+,k+)zk=Aij;Aij=Ai+1j;Ai+1j=zk;t1=bi;bi=bi+1;bi+1=t1;t2=Ai+1i;for(j=i;j<i+3;j+)Ai+1j=Ai+1j-Aij*t2/Aii;bi+1=bi+1-

37、bi*t2/Aii; x8=b8/A88;for(i=7;i>=0;i-)for(j=i+1;j<9;j+)t3=t3+Aij*xj;xi=(bi-t3)/Aii;t3=0;5.3運行結(jié)果*' D=VC + + EXERC5 + ewe"The JEes u 1of K Ai*e =1-1-G757992=2 27S7443=-2 合F95：L占"1=2 2930775-2-112&346=-6-4238377-1-3579238=0-6342449=-0-587266i*ess a.ny ko_y to continue5.4 MATLAB上機程序un cti on x=mgauss2(A,b,flag)if narginv 3,flag=0;e ndn=len gth(b);for k=1:( n-1)ap,p=max(abs(A(k: n, k);p=p+k-1;if p>kA(k p,:)=A(p k,:);b(k p,:)=b(p k,:);endm=A(k+1: n,k)/A(k,k);A(k+1: n,k+1: n)=A(k