1、-. z.*工學院數值分析考試基于Matlab的方法綜合應用報告班級：金融1121 *：姚婷婷*： 1124104129 成績：數理學院 2014年6月7日數值分析課程綜述報告前言：數值分析也稱計算方法，它與計算工具的開展密切相關。數值分析是一門為科學計算提供必需的理論根底和有效、實用方法的數學課程，它的任務是研究求解各類數學問題的數值方法和有關的理論。正文：第一章近似計算與誤差分析1、產生誤差的原因：模型誤差；觀測誤差；截斷誤差；舍入誤差。2、四則運算的誤差：加減法運算乘法運算除法運算：3、科學表示法、有效數字、近似值的精度任何一個實數都可以表示成如下的形式：其中：是正整數，是整數，如果是數

2、的近似值并且則稱該近似值具有位有效數字significant digit。此時，該近似值的相對誤差為另一方面，假設則，即：至少有位有效數字。例：取其近似值如下：*=3.14* =3.14159*=3.1415*=3.141第二章線性方程組在科學計算中，問題的本身就是求解線性方程組，許多問題的求解需要最后也歸結為線性方程組的求解，所以線性方程組的求解是科學計算中最常見的問題。對于線性方程組的求解一般有兩種方法：1) 直接法：高斯消去法；2) 間接法：各種迭代法。高斯消去法：求解思路：先消元，即按一定的規律逐步消去未知量，將方程組化為等價的上或下三角形方程組；然后進展回代，即由上三角形方程組逐個求

3、出；高斯列、全主元素消去法,及在消元的每一步選取列主元素列中絕對值最大的元取做主元素，計算步驟：消元過程：按列選主元、行交換、消元計算；回代過程；高斯列主元素消去法的MATLAB 實現：。第三章解線性方程組的迭代法通常逆矩陣不易求得，特別是對于大型的線性方程組，需要用迭代法求解。用迭代法求解線性方程組，要把線性方程組寫成等價的形式，右邊寫為迭代格式，如：2、關于迭代法收斂性的兩個重要結論：充分必要條件是：矩陣的譜半徑；充分條件是：矩陣的*個算子*數。3、線性方程組的迭代法主要有Jacobian迭代法，Gauss-Seidel迭代法。Jacobian 迭代法:Gauss-Seidel 迭代法:

4、(3.7)Jacobian 迭代法與G-S迭代法比擬：3.8式3.7和 (3.8) 說明：Gauss-Seidel 迭代法在計算第次迭代的第個分量時，及時地利用了在此步迭代中得到的新的迭代值：，由于第步的迭代值通常比第步的迭代值更接近方程組的準確解，所以，在Jacobian迭代法和GS迭代法都收斂的情況下，Gauss-Seidel迭代法的收斂速度比Jacobian迭代法的收斂速度高。例題：用MATLAB函數normrdn生成5階矩陣和向量分別構造線性方程組的Jacobi迭代格式和G-S迭代格式，并判斷收斂性。Jacobian迭代法和GS迭代法程序如下：clc;clear all;%1MbM=n

5、ormrnd(1,2,5)b=normrnd(1,2,5,1)%JacobianM1=D(L+U)f1=Dbrho=ma*(abs(eig(M1);R=1e-08; %switch sign(1-rho) case -1 disp(the Jocobian method is not applicable) otherwise *(:,1)=normrnd(0,9,5,1); k=1 while k=R k=k+1; else *=*(:,k+1); disp(Jacobian) IterN=k %Jacobian break end endend%Causs-Seidel M2=(D-L)U

6、f2=(D-L)brho=ma*(abs(eig(M2);R=1e-08;switch sign(1-rho) case -1 disp(the auss-seidel method is not applicable) otherwise *(:,1)=normrnd(0,9,5,1); k=1 while k=R k=k+1; else *=*(:,k+1) disp(Causs-Seidel) IterN=k break end endend第四章非線性代數方程組的數值解法：一、二分法：首先要確定適當的包含根的區間，這可以依據閉區間上連續函數的介值定理來確定，例如該方程：對于該方程所以該

7、方程至少有一個實根位于區間，圖像說明該區間中只含有一個實根；用表示方程在區間上的準確解，對于給定的精度要求，取區間的中點，并按下式進展判斷：1以為例，如果,沒有到達精度要求，令，并重復1的迭代過程；如果，則，必有，因為。即區間內的任何一點都可以作為方程的近根，特別地，可取做為近似解。二、牛頓迭代法：任取初始值上過點的切線方程為：與軸交于點，過點的切線方程為與軸交于點，如此下去得牛頓迭代公式：例題：考慮如下三階非線性方程組：其中取適當的迭代初值，用Newton迭代法求方程組的數值解程序：%Newton*=sym(*,clear);y=sym(y,clear);syms z; F=*2+y2+a2

8、*z2/2-a2; *+y-a; (2*-a)2+(2*y-a)2+a2*(z-b)/4; F*=diff(F,*,1);Fy=diff(F,y,1);Fz=diff(F,z,1); DF=F* Fy Fz; F=(*,y,z)*2+y2+a2*z2/2-a2;*+y-a; (2*-a)2+(2*y-a)2+a2*(z-b)/4; %NewtonFr=10-10; Er=10-10; *0=-a/4;-a;-a/2; *0=-a/4;-a;a/2;k=1; *(:,1)=*0;while norm(F(*(1,k),*(2,k),*(3,k)-0;0;0,2)=Fr tic; f(:,k)=F

9、(*(1,k),*(2,k),*(3,k); J=subs(DF,*,*(1,k); J=subs(J,y,*(2,k);J=subs(J,z,*(3,k); *(:,k+1)=*(:,k)-Jf(:,k); t(k)=toc; if norm(*(:,k+1)-*(:,k),2)Er breakend k=k+1;enddisp(Newton);disp(*(:,end);運行結果：Newton迭代法結果為：3.4291 0.46580.6535 第五章插值一、插值：插值是離散函數逼近的重要方法，利用它可通過函數在有限個點處的取值狀況，估算出函數在其他點處的近似值。二、常用差值法：拉格朗日L

10、agrange插值法、牛頓Newton插值法1、拉格朗日Lagrange插值法：拉格朗日插值多項式 = 1 * GB2 簡單的證明因為拉格朗日插值多項式的基函數有如下的性質：所以拉格朗日插值多項式滿足插值的條件。插值多項式 = 3 * GB2 拉格朗日插值法的缺乏在實際問題中，觀測的數據可能會不斷增加，如果用拉格朗日插值公式構造插值多項式，則，每當增加數據就要重新計算多項式的系數，由此增加許多不必要的計算工作量。2、三次樣條Spline插值 = 1 * GB2 插值條件要求要求所求的插值多項式三次樣條函數 = 1 * GB3 在每個區間，是次數不超過三次的多項式； = 2 * GB3 ，； =

11、 3 * GB3 在區間.上具有二階連續導數。例題：在夏季，較大湖泊的水體按深度被躍變層分為上部的變溫層和下部的均溫層。水體的分層化對環境工程中污染問題的研究具有重要的意義，例如，有機物的分解將導致被躍變層隔離的底部水體中氧急劇減少。按溫度隨水深的變化曲線，躍變層位于水深處：現有美國普拉特湖Platte Lake的一組數據：深度(): 02.34.99.113.718.322.927.2溫度(): 22.822.822.820.613.911.711.111.11. 試利用Lagrange差值方法求溫度隨水深近似變化函數表達式；2. 試利用三次樣條差值方法應用Matlab函數csape求溫度隨

12、水深近似變化函數表達式；3. 畫出插值函數和曲線，并與原始插值數據圖像作比擬程序代碼：函數文件程序：function Ln = my_Fun(*, *I, YI)if isa(*, sym) = 1; n = length(*I) - 1; Ln = 0; Pn = sym(ones(n + 1, 1);for k = 1 : n + 1for i = 1 : n + 1if i = k Pn(k) = Pn(k)*(* - *I(i)/(*I(k) - *I(i);else Pn(k) = Pn(k);endend Ln = Ln + Pn(k)*YI(k);endelse n = leng

13、th(*I) - 1; L = ones(n + 1, length(*); Ln = zeros(size(*);for k = 1 : n + 1for i = 1 : n + 1if i = k L(k,:) = L(k,:).*(* - *I(i)./(*I(k) - *I(i);else L(k,:) = L(k,:);endend Ln = Ln + L(k, :).*YI(k);endEnd主文件程序：clc;clear all;close all;%1.Lagranget = sym(t, real);*=0 2.3 4.9 9.1 13.7 18.3 22.9 27.2;T=

14、22.8 22.8 22.8 20.6 13.9 11.7 11.1 11.1;*=linspace(0,27.2,275);Ln=my_Fun(*,*,T);figure(1)plot(*,T,r*,markersize,15) *label(*);ylabel(T);title();pause(3)hold onplot(*,Ln,b-,linewidth,3);*label(*);ylabel(T);title(Lagrange);%I_poly = my_Fun(t, *, T);I_poly = simple(I_poly)I_poly = sym2poly(I_poly)%2. f

15、igure(2)hold ony=interp1(*,T,*,spline);plot(*,y,b-,linewidth,3); %*label(*);ylabel(T);title();PP=csape(*,T,2,2,0,0); Coefs = PP.coefs%3.figure(3)plot(*,Ln,-,color,1, 0, 0,LineWidth, 3)*label(*);ylabel(T);title(Lagrange);hold onpause(3)fnplt(PP,b,3,0,28) %*label(*);ylabel(T);title();pause(3)plot(*,T,

16、-,linewidth,3,markersize,10)運行結果:I_poly =(0*t7 - 00*t6 + 0000*t5 - 35000*t4 + *t3 - *t2 + *t + 28)/0I_poly = Columns 1 through 6 0.0000 -0.0000 0.0013 -0.0193 0.1268 -0.3685 Columns 7 through 8 0.3770 22.8000Coefs = 0.0022 -0.0000 -0.0115 22.8000 -0.0092 0.0150 0.0230 22.8000 -0.0114 -0.0565 -0.0850

17、 22.8000 0.0297 -0.2004 -1.1640 20.6000 -0.0153 0.2099 -1.1200 13.9000 0.0017 -0.0014 -0.1605 11.7000 -0.0017 0.0223 -0.0640 11.1000圖示如下：第六章最小二乘擬合與最正確逼近一、最小二乘擬合加權最小二乘法逼近準則：最小二乘逼近多項式必須滿足如下必要條件：即滿足法方程組:例題：下面是一處地質巖層斷面高程深度的測量數據。水平距離():00.201.002.103.505.006.807.509.0011.212.0高度():h 1.641.581.681.841.580

18、.860.390.310.390.770.86試利用最小二乘法求滿足即誤差不超過的最低次數的擬合多項式,寫出該多項式的表達式;畫出數據散點圖和該多項式曲線.程序：clc;clear all;close all;*=linspace(0,12,12);*=0 200 1000 2100 3500 5000 6800 7500 9000 11200 12000;h=1.64 1.58 1.68 1.84 1.58 0.86 0.39 0.31 0.39 0.77 0.86;plot(*,h,*,markersize,8) %畫出給出數據的散點圖*label(*);ylabel(h);title(*

19、h);figure(1)%1. n=3;P, S = polyfit(*, h, n)Pm = polyval(P, *)R = sqrt(sum(h - Pm).2) %t = linspace(*(1), *(end), 12);Poly = polyval(P, t);figure(2)plot(*, h, ro, t, Poly, LineWidth, 3, markersize, 8)set(gca,FontSize,18)legend(The Data, The Fitting Curve, 1)title(Curve Fitting by Least Square Appro*i

20、mation, fontsize, 18)第七章微積分的數值方法一、數值微分如果給定函數的關系式比擬復雜或者未知，而僅僅知道在個相異點，處的函數值，則，我們可以利用函數的插值多項式的導數作為函數導數的近似值因而有這里需要說明一點的是，盡管和的函數值可能相差不多，但是它們的導數有可能相差很大。二、數值積分數值積分方法最大的優點是將復雜的函數積分轉化為相對簡單的加、減、乘、除運算。對于給定的函數，如果的函數關系式比擬復雜或.未知，而僅僅知道該函數在個點處的函數值，則可以利用函數的插值多項式的積分作為函數在上的定積分的近似值。1、Newton-Cotes求積公式根據Lagrange插值多項式有令得N

21、ewton-Cotes求積公式特別地，當取插值節點為時有 = 1 * GB3 兩點公式梯形公式： = 2 * GB3 三點公式Simpson公式：例題：下面是一處地質巖層斷面上部邊緣的深度測量數據。水平距離():00.201.002.103.505.006.807.509.0011.212.0深度(): 1.641.581.681.841.580.860.390.310.390.770.86表1 = 1 * GB2 試利用復化的梯形求積法求該組數據所在曲線與基準線軸在*圍內所圍成圖形面積.畫出數據散點圖和圖形的示意圖. = 2 * GB2 試利用復化的Simpson求積法求該組數據所在曲線與基

22、準線軸在*圍內所圍成圖形面積.畫出數據散點圖和圖形的示意圖.解答如下：程序：% *I = 0 0.20 1.00 2.10 3.50 5.00 6.80 7.50 9.00 11.2 12.00;YI = -1.64 1.58 1.68 1.84 1.58 0.86 0.39 0.31 0.39 0.77 0.86;n = length(*I) - 1; %n=10M = 60;* = linspace(*I(1), *I(end), M +1); % 1.1 T1 = 0; %for k = 2 : n + 1 T1 = T1 + (YI(k) + YI(k-1)*(*I(k) - *I(k

23、-1)/2;enddisp()T_Trape = T1 %figure(1) %set(gca,FontSize, 20)patch(*I(1), *I, *I(end), 0, YI, 0, 0.5 1 0.5) %hold onplot(*I, YI, o-, *I, 0*I, k, LineWidth, 3, markersize, 8)title()*label(*)ylabel(y)運行結果：復化的梯形求積法求得的圖形面積：T_Trape = -11.6090圖像：程序：函數M文件：function Ln = Lagrange_Fun_01(*, *I, YI)if isa(*, s

24、ym) = 1; n = length(*I) - 1; Ln = 0; Pn = sym(ones(n + 1, 1);for k = 1 : n + 1for i = 1 : n + 1if i = k Pn(k) = Pn(k)*(* - *I(i)/(*I(k) - *I(i);else Pn(k) = Pn(k);endend Ln = Ln + Pn(k)*YI(k);endelse n = length(*I) - 1; L = ones(n + 1, length(*); Ln = zeros(size(*);for k = 1 : n + 1for i = 1 : n + 1

25、if i = k L(k,:) = L(k,:).*(* - *I(i)./(*I(k) - *I(i);else L(k,:) = L(k,:);endend Ln = Ln + L(k, :).*YI(k);endend主程序：% 1.2 1.3%32S2 = 0;%Simpsonfor k = 1: 2: n - 1 S2 = S2 + (*I(k+2) - *I(k)*(YI(k) + 4*YI(k+1) + YI(k+2)/6; *(k, :) = linspace(*I(k), *I(k+2), 20); Y(k, :) = Lagrange_Fun_01(*(k,:), *I(k

26、 : k+2), YI(k : k+2);enddisp(Simpson)S_Simpson = S2 %Simpsonfigure(2)set(gca,FontSize, 20)hold onfor k = 1: 2 : n-1 patch(*I(k), *(k,:), *I(k +2), 0, Y(k, :), 0, 0.5 1 0.5) plot(*(k, :), Y(k, :), LineWidth, 3)endplot(*I, YI, o, *I, 0*I, k, LineWidth, 3, markersize, 8)title(Simpson)*label(*)ylabel(y)

27、運行結果：Simpson復化求積法求得的圖形面積：S_Simpson = -11.9128圖像：第八章微分方程組初值問題數值方法引言：微分方程組是科學研究和工程應用中最常用的數學模型之一。但是基絕大多數的常微分方程和常微分方程組得不到解析解和實際應用中往往只需要知道常微分方程組的解在*些點處的函數值這兩點原因，絕大多數在實際中遇到的常微分方程和常微分方程組得不到解析解，我們可以采用下面將介紹的常微分方程組的初值問題的數值解法，就可以到達這一目的。Euler 方法最簡單的數值解法1、Euler 法假設要求在點,，處初值問題4的解的近似值。首先對式4的兩端積分，得對于該式的右邊，如果被積函數用積分

28、下限處的函數值代替被積函數作積分，則有進而得到下式給出的遞推算法Euler 方法2、Euler 方法的改良改良的Euler 方法Euler 方法的精度高，其原因在于：在推導Euler 方法時，我們是用待求解函數在一點處的變化率1代替在區間上的平均變化率；而在推導改良的Euler 方法時，我們是用待求解函數在兩點處變化率的平均值2代替在區間上的平均變化率；顯然，通常式2比式1更接近于在區間上的平均變化率。由此，人們受到啟發：適當地選取區間上函數多個點處的變化率，用它們加權平均值代替在區間上的平均變化率，近似解的精度應更高。RungeKutta法： = 1 * GB3 二階的RungeKutta法 = 2 * GB3 三階的RungeKutta法 = 3 * GB3 四階的RungeKutta法RungeKutta法是按選取區間上函數變化率數目的多少和截斷誤差的階數來區分的一系列方法。例題：給定如下常微分方程Cauchy 問題：分別取參數，初值和步長：，。取，在區間上分別利用Euler 法，4階R-K公式和matlab的ode45求問題1的數值解；程序：clcclear alla=1;b=1;F = (t, y) -a*y-t2*e*p(-0.5*t)*sin(t)+b;%matlabode45y0 = -1;Tf = 2*pi;Tspan = 0,