數(shù)值計算實驗

1、重 慶 大 學(xué)學(xué) 生 實 驗 報 告實驗課程名稱 數(shù)值計算 開課實驗室 數(shù)學(xué)與統(tǒng)計學(xué)院 學(xué) 院 數(shù)學(xué)與統(tǒng)計學(xué)院 年級 X級專業(yè)班 數(shù)學(xué)與應(yīng)用數(shù)學(xué)一班 學(xué) 生 姓 名 X 學(xué) 號 X 開 課 時 間 2012 至 2013 學(xué)年第 一 學(xué)期總 成 績教師簽名課程名稱數(shù)值計算實驗項目名 稱Gauss消元法實驗項目類型驗證演示綜合設(shè)計其他指導(dǎo)教師成 績一、實驗?zāi)康模?）高斯列主元消去法求解線性方程組的過程（2）熟悉用迭代法求解線性方程組的過程（3）設(shè)計出相應(yīng)的算法，編制相應(yīng)的函數(shù)子程序二、實驗內(nèi)容分別用高斯列主元消去法 ,Jacobi迭代法,Gauss-Saidel迭代法,超松弛迭代法求解線性方程組

2、三、實驗原理（1）高斯列主元消去法：設(shè)有方程組AX=b，將矩陣的初等行變換作用于方程組大的增廣矩陣B=A:b，將A變換成一個上三角矩陣，然后對方程組進行求解。（2）Jacobi迭代法：設(shè)方程組Ax=b 滿足aii 0, 將方程組變形為: x=Bx+f, 則雅可比(Jacobi)迭代法是指x(k+1)=Bx(k)+f，即由初始解逐步迭代即可得到方程組的解。（3）Gauss-Saidel迭代法: 在雅可比迭代法計算過程中,在迭代的每一步計算過程中是用的全部分量來計算的所有分量,顯然在計算第i個分量時,已經(jīng)計算出的最新分量沒有被利用，對這些最新計算出來的第次近似的分量加以利用,就得到所謂解方程組的G

3、auss-Seidel迭代法。（4）超松弛迭代法：在Gauss-Saidel迭代法的基礎(chǔ)上，對每次迭代是的改變量乘以一個因子w，是收斂速度加快。四、實驗步驟（1）高斯列主元消去法：根據(jù)實驗原理編寫matlab程序，首先判斷解的個數(shù)，再分情況進行列主元消去求解。Matlab程序如下：clcA=input('A');b=input('b');B=A b' n=length(b);RA=rank(A);RB=rank(B);c=RB-RA;if c>0, disp('方程組無解') returnendif RA=RB if RA=n fp

4、rintf('方程組有唯一解',n) X=zeros(n,1); for p=1:n-1 t=find(abs(B(p:end,p)=max(abs(B(p:end,p)+p-1; if abs(B(t,p)=abs(B(p,p) l=B(t,:); B(t,:)=B(p,:); B(p,:)=l; end for k=p+1:n m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1); end end b=B(1:n,n+1); A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1

5、 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); end else disp('方程組有無窮多解') return endend disp('X');X(2) Jacobi迭代法: 根據(jù)實驗原理編寫matlab程序，調(diào)用Jacobi_iterative函數(shù)，首先給定誤差限，確定最大迭代次數(shù)，確認(rèn)矩陣A的階，然后進行迭代過程，最后輸出結(jié)果。Matlab程序如下：A=2 10 0 -3;-3 -4 -12 13;1 2 3 -4;4 14 9 -13;b=10 5 -2 7'e=0.000001; %控制誤差n=max(

6、size(A); %測定維數(shù)for i=1:n if A(i,i)=0 '對角元為零，不能求解' return endendx=zeros(n,1) %設(shè)置初始解k=0; %預(yù)設(shè)迭代次數(shù)為0kend=50 %最大迭代次數(shù)為50r=1; %前后項之差的無窮范數(shù)，初始值設(shè)為1 while k<=kend & r>e %達到預(yù)定精度或迭代超過50次退出計算 x0=x; %記下前次近似解 for i=1:n s=0; for j=1:i-1 s=s+A(i,j)*x0(j); end for j=i+1:n s=s+A(i,j)*x0(j); end x(i)=b(

7、i)/A(i,i)-s/A(i,i); end r=norm(x-x0,inf); %重新計算前后項之差的無窮范數(shù) k=k+1;endif k>kend '迭代不收斂，失敗 'else '求解成功' x kend（3）Gauss-Saidel迭代法：根據(jù)實驗原理編寫matlab程序，調(diào)用Gauss_Seidel_iterative函數(shù)，首先給定誤差限，確定最大迭代次數(shù)，確認(rèn)矩陣A的階，然后進行迭代過程，最后輸出結(jié)果。Matlab程序如下：function x=Gauss_Seidel_iterative(A,b) x0=zeros(1,length(b);

8、 tol=10(-3); N=10; n,n=size(A); k=1; while k<=N x(1)=(b(1)-A(1,2:n)*x0(2:n)')/A(1,1); for i=2:n x(i)=(b(i)-A(i,1:i-1)*x(1:i-1)'-A(i,i+1:n)*x0(i+1:n)')./A(i,i); end if max(abs(x-x0)<=tol fid = fopen('G_S_iter_result.txt', 'wt'); fprintf(fid,'方程組的解'); fprintf(

9、fid,'迭代次數(shù): %d次',k); fprintf(fid,'x的值'); fprintf(fid, '%12.8f', x); break; end k=k+1; x0=x; end if k=N+1 fid = fopen('G_S_iter_result.txt', 'wt'); fprintf(fid,'方程組的解'); fprintf(fid,'迭代次數(shù): %d次',k); fprintf(fid,'超過最大迭代次數(shù)'); fclose(fid); en

10、d（4）超松弛迭代法：根據(jù)實驗原理編寫matlab程序，調(diào)用SOR_iterative(A,b)函數(shù)，首先給定誤差限，確定最大迭代次數(shù)，確認(rèn)矩陣A的階，給定松弛因子，然后進行迭代過程，最后輸出結(jié)果。Matlab程序如下：function x=SOR_iterative(A,b) x0=zeros(1,length(b); tol=10(-2); N=1000; n,n=size(A); w=1; k=1; while k<=N x(1)=(b(1)-A(1,2:n)*x0(2:n)')/A(1,1); for i=2:n x(i)=(1-w)*x0(i)+w*(b(i)-A(i,

11、1:i-1)*x(1:i-1)'-A(i,i+1:n)*x0(i+1:n)')/A(i,i); end if max(abs(x-x0)<=tol fid = fopen('SOR_iter_result.txt', 'wt'); fprintf(fid,'方程組的解'); fprintf(fid,'迭代次數(shù): %d次',k); fprintf(fid,'x的值'); fprintf(fid, '%12.8f', x); break; end k=k+1; x0=x; end

12、if k=N+1 fid = fopen('SOR_iter_result.txt', 'wt'); fprintf(fid,'方程組的解'); fprintf(fid,'迭代次數(shù): %d次',k); fprintf(fid,'超過最大迭代次數(shù)'); fclose(fid);end五、 實驗結(jié)果1高斯列主元消元法：A=2 10 0 -3；-3 -4 -12 13；1 2 3 -4；4 14 9 -13；b=10 5 -2 7'；計算結(jié)果：>> GaussX= 0.999999999999999

13、2.000000000000000 2.999999999999999 3.999999999999999.2. Jacobi迭代法：A=2 10 0 -3;-3 -4 -12 13;1 2 3 -4; 4 14 9 -13；b=10 5 -2 7A = 2 10 0 -3 -3 -4 -12 13 1 2 3 -4 4 14 9 -13b = 10 5 -2 7>> Jacobi_iterative(A,b)ans = 1.0e+58 * -2.7656 -1.9132 -0.8584 0.88953.Gauss-Saidel迭代法：>> A=2 10 0 -3;-3

14、 -4 -12 13;1 2 3 -4; 4 14 9 -13,b=10 5 -2 7A = 2 10 0 -3 -3 -4 -12 13 1 2 3 -4 4 14 9 -13b = 10 5 -2 7>> Gauss_Seidel_iterative(A,b)ans = NaN NaN NaN NaN（4）超松弛迭代法：>> A=2 10 0 -3;-3 -4 -12 13;1 2 3 -4; 4 14 9 -13,b=10 5 -2 7A = 2 10 0 -3 -3 -4 -12 13 1 2 3 -4 4 14 9 -13b = 10 5 -2 7>&g

15、t; SOR_iterative(A,b)ans = NaN NaN NaN NaN六、實驗中遇到的問題及解決方法1）實驗中需要注意許多細節(jié)，如松弛因子的取值區(qū)間應(yīng)該是（0,2），最大迭代次數(shù)不能太大，以及高斯列主元消去法中方程組解的個數(shù)應(yīng)在消元前予以確認(rèn)。2）一開始寫程序的時候沒有加上修正操作行值的語句，導(dǎo)致了錯誤結(jié)果的輸出，通過反復(fù)檢查，加上了語句“m=m+k-1”來對操作行的值進行修正，得出了正確的結(jié)果。七、實驗結(jié)論通過實驗，我發(fā)現(xiàn)高斯列主元消去法可以得到方程組的解為X=1 2 3 4,其他兩種種方法求出方程組的解都是不可行的。而且可以得出Gauss消元法解線性方程組是適合在計算機上自動

16、計算的一種標(biāo)準(zhǔn)化方法。教師簽名年 月 日課程名稱數(shù)值計算實驗項目名 稱插值方法實驗項目類型驗證演示綜合設(shè)計其他指導(dǎo)教師何光輝成 績一、實驗?zāi)康模?1)學(xué)會拉格朗日插值、牛頓插值等基本方法(2)設(shè)計出相應(yīng)的算法，編制相應(yīng)的函數(shù)子程序(3)會用這些函數(shù)解決實際問題二、實驗內(nèi)容：（1）設(shè)計拉格朗日插值算法，編制并調(diào)試相應(yīng)的函數(shù)子程序（2）設(shè)計牛頓插值算法，編制并調(diào)試相應(yīng)的函數(shù)子程序（3）給定函數(shù)四個點的數(shù)據(jù)如下：X1.12.33.95.1Y3.8874.2764.6512.117試用拉格朗日插值確定函數(shù)在x=2.101，4.234處的函數(shù)值。（4）已知用牛頓插值公式求的近似值。三、實驗原理：（1）拉

17、格朗日插值:經(jīng)過n+1個點（x0,y0）(x1,y1).(xn,yn)構(gòu)造一個n次多項式，形如,使得 成立，其中為差值基函數(shù)。(2)牛頓差值：四、實驗步驟（1）拉格朗日插值:根據(jù)實驗原理編寫matlab程序，通過表格給出的已知數(shù)據(jù)得出差值多項式，再代入定點求得函數(shù)值。Matlab程序如下：function y=lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(

18、k)+s; end y(i)=s;end（2）牛頓差值:根據(jù)實驗原理編寫matlab程序，通過已知數(shù)據(jù)計算差商矩陣求得牛頓差值多項式，再代入定點求出函數(shù)值。Matlab程序如下：function newtfunc=newton(x,y)n=length(x);d=zeros(n,n);d(:,1)=y'syms s;for j=2:n for i=j:n d(i,j)=(d(i,j-1)-d(i-1,j-1)/(x(i)-x(i-j+1); endendtemp=1;newtfunc=d(1,1);for i=2:n temp=temp*(t-x(i-1); newtfunc=newt

19、func+temp*d(i,i);endunction f=f(t)f=t/3 - (t - 1)*(t - 4)/60 + 2/3;end五、 實驗結(jié)果:（1）>> x0=1.1 2.3 3.9 5.1,y0=3.887 4.276 4.651 2.117,x=2.101x0 = 1.1000 2.3000 3.9000 5.1000y0 = 3.8870 4.2760 4.6510 2.1170x = 2.1010>> lagrange(x0,y0,x)ans = 4.1457>> x=4.234x = 4.2340>> lagrange(x

20、0,y0,x)ans =4.3007即函數(shù)在x=2.101，4.234處的函數(shù)值分別為4.2340 4.3007（2）>> x=1 4 9,y=1 2 3x = 1 4 9y = 1 2 3>> newton(x,y) ans =t/3 - (t - 1)*(t - 4)/60 + 2/3 >> f(t)ans =2.2667即的近似值為2.2667六、實驗中遇到的問題及解決辦法實驗中牛頓差值多項式無法直接得出函數(shù)值，解決辦法是通過定義一個新函數(shù)將定點代入多項式中求出函數(shù)值，算法還可以優(yōu)化，使得牛頓差值多項式的值直接得出。七、 實驗結(jié)論通過拉格朗日差值和牛頓

21、差值得出了多項式在定點的值，這兩種算法的matlab程序調(diào)試正常，可以用這兩種算法解決實際問題。教師簽名年 月 日課程名稱數(shù)值計算實驗項目名 稱數(shù)值微積分實驗項目類型驗證演示綜合設(shè)計其他指導(dǎo)教師何光輝成 績一、實驗?zāi)康模海?）學(xué)會復(fù)化梯形、復(fù)化辛浦生求積公式的應(yīng)用（2）設(shè)計出相應(yīng)的算法，編制相應(yīng)的函數(shù)子程序（3）會用這些函數(shù)解決實際問題二、實驗內(nèi)容：（1）設(shè)計復(fù)化梯形公式求積算法，編制并調(diào)試相應(yīng)的函數(shù)子程序（2）設(shè)計復(fù)化辛浦生求積算法，編制并調(diào)試相應(yīng)的函數(shù)子程序（3）分別用復(fù)化梯形公式和復(fù)化辛浦生公式計算定積分取n=2，4，8，16，精確解為0.9460831三、實驗原理：（1）復(fù)化梯形：每個

22、小區(qū)間xk，xk+1上梯形公式是：整個區(qū)間上的復(fù)化梯形公式是： 即：（2）復(fù)化辛浦生：記每個小區(qū)間xk，xk+1的中點為xk+1/2，該小區(qū)間上的辛甫生公式為：整個區(qū)間上的復(fù)化辛甫生公式是：即：四、實驗步驟：（1）復(fù)化梯形：根據(jù)實驗原理編寫matlab程序，首先確定函數(shù)f，通過給定區(qū)間及n計算步長，再計算端點的函數(shù)值，再進行累加得出積分值。Matlab程序：function t=tn(f,a,b,n)f=inline(f);h=(b-a)/n;syms x t=(limit(f(x),a)+limit(f(x),b)*h/2; for k=1:n-1 t=t+h*f(a+k*h);endt=v

23、pa(t);（2）復(fù)化辛浦生：根據(jù)實驗原理編寫matlab程序，程序與復(fù)化梯形的程序相似，只是最后的累加公式略有不同。Matlab程序：function t=tn(f,a,b,n)f=inline(f);h=(b-a)/n;syms x t=(limit(f(x),a)+limit(f(x),b)*h/6; for k=1:n-1 t=t+h*f(a+k*h)/3+2*h*f(a+(2*k-1)*h/2)/3;endt=t+2*(b-a)/n*f(a+(2*n-1)*(b-a)/n/2)/3t=vpa(t);五、實驗結(jié)果：分別取n=2，4，8，16（1）>> tn('sin

24、(x)/x',0,1,2) t = sin(1)/12 + 94679825489681579/108086391056891904 ans = 0.94608693395179367210035361043684 >> tn('sin(x)/x',0,1,4) t = sin(1)/24 + 196938164598932969/216172782113783808 ans = 0.94608331088847189475438310411692 >> tn('sin(x)/x',0,1,8) t = sin(1)/48 + 4

25、01455528529665361/432345564227567616 ans = 0.94608308538494763586154759765045 >> tn('sin(x)/x',0,1,16) t = sin(1)/96 + 202622585428064699/216172782113783808 ans = 0.94608307130556208489984050743676（2）>> tn('sin(x)/x',0,1,2) t = sin(1)/12 + 94679825489681579/10808639105689

26、1904 ans = 0.94608693395179367210035361043684 >> tn('sin(x)/x',0,1,4) t = sin(1)/24 + 196938164598932969/216172782113783808 ans = 0.94608331088847189475438310411692 >> tn('sin(x)/x',0,1,8) t = sin(1)/48 + 401455528529665361/432345564227567616 ans = 0.946083085384947635861

27、54759765045 >> tn('sin(x)/x',0,1,16) t = sin(1)/96 + 202622585428064699/216172782113783808 ans = 0.94608307130556208489984050743676六、實驗中遇到的問題及解決辦法被積函數(shù)在端點沒有定義，造成計算發(fā)生錯誤，解決方法是通過取端點處的函數(shù)極限值來代替函數(shù)在端點的值。七、實驗結(jié)論兩種算法的matlab程序調(diào)試正常，可以用來解決實際問題，通過不同的算法都得出了積分值，而且兩種算法得出的積分值與精確值很接近，所以，復(fù)化梯形公式和復(fù)化辛浦生公式的精度都

28、很高。 教師簽名年 月 日 課程名稱數(shù)值計算實驗項目名 稱常微分方程的數(shù)值解法實驗項目類型驗證演示綜合設(shè)計其他指導(dǎo)教師何光輝成 績一、 實驗?zāi)康模海?）學(xué)會四階龍格-庫塔方法的使用（2）設(shè)計出相應(yīng)的算法，編制相應(yīng)的函數(shù)子程序（3）會用這些函數(shù)解決實際問題二、 實驗內(nèi)容：（1）分別取h=0.05，N=10；h=0.025，N=20；h=0.01，N=50，用四階龍格-庫塔方法求解微分方程初值問題：y=-50y，y(0)=10（2）某跳傘者在t=0時刻從飛機上跳出，假設(shè)初始時刻的垂直速度為0，且跳傘者垂直下落。已知空氣阻力為F=cv2，其中c為常數(shù)，v為垂直速度，向下方方向為正。寫出此跳傘者的速度

29、滿足的微分方程；若此跳傘者的質(zhì)量為M=70kg，且已知c=0.27kg/m，利用四階龍格-庫塔公式計算t<=20s的速度（取h=0.1s）三、實驗原理：龍格-庫塔(Runge-Kutta)方法是一種在工程上應(yīng)用廣泛的高精度單步算法。四階經(jīng)典龍格庫塔算法的公式：yi+1=yi+h*( K1+ 2*K2 +2*K3+ K4)/6K1=f(xi,yi)K2=f(xi+h/2,yi+h*K1/2)K3=f(xi+h/2,yi+h*K2/2)K4=f(xi+h,yi+h*K3)四、 實驗步驟：（1）根據(jù)實驗原理編寫matlab程序，首先確定區(qū)間，給定步長，通過給定的初值計算并輸出函數(shù)值，其中，要使

30、用不同的步長進行計算。Matlab程序：function l=m(p)clc,clearx0=0;xn=0.5;y0=10;h=0.05;y,x=g(x0,xn,y0,h); n=length(x);fprintf(' i x(i) y(i)n');for i=1:nfprintf('%2d %4.4f %4.4fn',i,x(i),y(i);endfunction z=f(x,y)z=50*y;function y,x=g(x0,xn,y0,h) x=x0:h:xn;n=length(x);y1=x;y1(1)=y0;for i=1:n-1 K1=f(x(i)

31、,y1(i); K2=f(x(i)+h/2,y1(i)+h/2*K1); K3=f(x(i)+h/2,y1(i)+h/2*K2); K4=f(x(i)+h,y1(i)+h*K3); y1(i+1)=y1(i)+h/6*(K1+2*K2+2*K3+K4);end y=y1;（2）首先建立數(shù)學(xué)模型，通過物理知識可得，dv/dt =(Mg-cv2)/M,問題轉(zhuǎn)化為求解常微分方程，編寫相應(yīng)的matlab程序，采用ode45直接求解。Matlab程序：function f=fun(t,v)f=(9.8*70-0.27*v2)/70;endt,v=ode23('fun',0,20,1);p

32、lot(t,v,'o-');t,v=ode45('fun',0,20,1); plot(t,v,'*-')五、實驗結(jié)果：(1) h=0.05，N=10 i x(i) y(i) 1 0.0000 10.0000 2 0.0500 108.5677 3 0.1000 1178.6947 4 0.1500 12796.8186 5 0.2000 138932.1265 6 0.2500 1508354.2586 7 0.3000 16375856.5211 8 0.3500 177788921.4496 9 0.4000 1930213576.8840

33、10 0.4500 20955886463.618611 0.5000 227513256944.8590h=0.025，N=20i x(i) y(i) 1 0.0000 10.0000 2 0.0250 34.5850 3 0.0500 119.6120 4 0.0750 413.6775 5 0.1000 1430.7019 6 0.1250 4948.0770 7 0.1500 17112.9049 8 0.1750 59184.9148 9 0.2000 204690.796510 0.2250 707922.320111 0.2500 2448346.578912 0.2750 84

34、67597.079213 0.3000 29285151.422014 0.3250 101282581.797715 0.3500 350285413.512416 0.3750 1211460734.330117 0.4000 4189832217.412118 0.4250387719 0.4500 50115401135.438120 0.4750 173323919063.627121 0.5000 599440097034.9956h=0.01，N=50i x(i) y(i) 1 0.0000 10.0000 2 0.0100 16.4844 3 0.02

35、00 27.1735 4 0.0300 44.7938 5 0.0400 73.8397 6 0.0500 121.7201 7 0.0600 200.6480 8 0.0700 330.7557 9 0.0800 545.230210 0.0900 898.777911 0.1000 1481.579112 0.1100 2442.290613 0.1200 4025.963514 0.1300 6636.549115 0.1400 10939.936416 0.1500 18033.801517 0.1600 29727.594618 0.1700 49004.081819 0.1800

36、80780.166120 0.1900 133161.055021 0.2000 219507.676622 0.2100 361844.685723 0.2200 596478.349124 0.2300 983257.278525 0.2400 1620838.170126 0.2500 2671850.421027 0.2600 4404378.428428 0.2700 7260342.565529 0.2800 11968220.947830 0.2900 19728864.218731 0.3000 32521799.610532 0.3100 53610154.045433 0.

37、3200 88372988.309134 0.3300 145677347.915835 0.3400 240140003.205036 0.3500 395855786.533337 0.3600 652543523.113538 0.3700 1075677213.882339 0.3800 1773186657.259240 0.3900 2922987380.325641 0.4000 4818362009.755642 0.4100 7942768625.456443 0.4200025844 0.4300 21583252073.605145 0.4400

38、 35578642090.083446 0.4500 58649167820.371847 0.4600 96679487578.894248 0.4700 159370092805.833449 0.4800 262711637359.616050 0.4900 433063714709.992151 0.5000 713878467217.2526(2)六、 實驗中遇到的問題及解決辦法在（2）中求解V時，由于步長較小，區(qū)間較大，采用（1）的方法會得出很多個解，不易觀察，所以采用ode45直接輸出圖像。七、實驗結(jié)論四階龍格-庫塔方法可以接出常微分方程。而且四階龍格-庫塔方法可以很好地求解常微分方程中，并且能夠在較好地運用于解決實際問題，是一種經(jīng)典而實用的方法。 教師簽名年 月 日 課程名稱數(shù)值計算實驗項目名 稱估計水塔的水流量實驗項目類型驗證演示綜合設(shè)計其他指導(dǎo)教師何光輝成 績一、實驗?zāi)康模海?）學(xué)會對實際問題的分析方法（2）學(xué)會利用所學(xué)的知識解決實際問題（3）設(shè)計出相應(yīng)的算法，編制相應(yīng)的應(yīng)用程序二、實驗內(nèi)容：某居民區(qū)，其自來水是有一個圓柱形水塔提供，水塔高12.2m，塔的直徑為17.4m，水塔是由水泵根據(jù)水塔中的水位自動加水，一般水泵每天工作兩次。按照設(shè)計，當(dāng)水塔中的水位降低至最低水位，約8.2m時，水泵自動啟動加水。當(dāng)水位升至最高水位，約10.8m