數(shù)值分析課程設(shè)計(jì)

1、課程設(shè)計(jì)（論文）任務(wù)書學(xué)院理學(xué)院專業(yè)信息與計(jì)算科學(xué)學(xué)生姓名班級(jí)學(xué)號(hào)課程名稱數(shù)值分析課程設(shè)計(jì)課程設(shè)計(jì) （論文）題目數(shù)值分析算法設(shè)計(jì)與應(yīng)用案例設(shè)計(jì)要求（技術(shù)參數(shù)）：1、熟練掌握Matlab軟件。2、運(yùn)用Matlab軟件進(jìn)行各種數(shù)值算法編程。3、每個(gè)同學(xué)至少做四個(gè)算法編程實(shí)驗(yàn)，每種算法最少有一個(gè)實(shí)例。4、每種算法在關(guān)鍵語句中要有注釋。設(shè)計(jì)任務(wù)：實(shí)驗(yàn)一：LU分解法解線性方程組實(shí)驗(yàn)二：Lagrange插值法數(shù)值求解實(shí)驗(yàn)三：龍貝格求積公式求數(shù)值積分實(shí)驗(yàn)四：用Runge-Kutta方法求常微分方程數(shù)值解實(shí)驗(yàn)五：線性規(guī)劃任務(wù)分配問題實(shí)驗(yàn)六：求電阻問題計(jì)劃與進(jìn)度安排：第一周：查閱資料 學(xué)習(xí)使用Matlab軟件完

2、成實(shí)驗(yàn)一、實(shí)驗(yàn)二第二周：完成實(shí)驗(yàn)三、實(shí)驗(yàn)四、實(shí)驗(yàn)五 攥寫課程設(shè)計(jì)報(bào)告 答辯成績：指導(dǎo)教師（簽字）：年 月日專業(yè)負(fù)責(zé)人（簽字）：年 月日主管院長（簽字）年 月日目錄實(shí)驗(yàn)一 LU 分解法解線性方程組 11.1 實(shí)驗(yàn)?zāi)康呐c要求 11.2 實(shí)驗(yàn)基本原理 11.3 算法程序 21.4 實(shí)例及結(jié)果分析 3實(shí)驗(yàn)二 Lagrange 插值法數(shù)值求解 42.1 實(shí)驗(yàn)?zāi)康呐c要求 42.2 實(shí)驗(yàn)基本原理 42.3 算法程序 42.4 實(shí)例及結(jié)果分析 5實(shí)驗(yàn)三 龍貝格求積公式求數(shù)值積分 53.1 實(shí)驗(yàn)?zāi)康呐c要求 53.2 實(shí)驗(yàn)基本原理 53.3 算法程序 63.4 實(shí)例及結(jié)果分析 7實(shí)驗(yàn)四 用 Runge-Kutta

3、方法求常微分方程數(shù)值解 74.1 實(shí)驗(yàn)?zāi)康呐c要求 74.2 實(shí)驗(yàn)基本原理 74.3 算法程序 84.4 實(shí)例及結(jié)果分析 8實(shí)驗(yàn)五 線性規(guī)劃任務(wù)分配問題 115.1 問題提出 115.2 模型分析 115.3 模型建立 125.4 模型求解 12實(shí)驗(yàn)六 求電阻問題 136.1 問題提出 136.2 問題分析與模型建立 135.3 模型求解 146.4 結(jié)果分析 15j -1Uij 二 aji - l jkUkik4i A(2) yi 二 bi - lik ykk nXi 二(yi -、 Uik Xk )/Uii實(shí)驗(yàn)一 LU分解法解線性方程組1.1實(shí)驗(yàn)?zāi)康呐c要求1了解LU分解法解線性方程組的基本原

4、理；2. 熟悉計(jì)算方法的技巧和過程，能用 LU分解法解實(shí)際問題;3. 用matlab實(shí)現(xiàn)LU分解。12實(shí)驗(yàn)基本原理1若一個(gè)線性方程組系數(shù)矩陣為n階方陣A且各階順序主子式均不為0則A 的LU分解存在且唯一。2. 在滿足1的條件下課推導(dǎo)得出以下公式j(luò)(1) lij 二佝 - l ikUij ) / UijkA3. 公式（1）用于求解矩陣L、U，公式（2）用于會(huì)帶求解y、x。從公式中 可以看出：L對(duì)角線上元素為1， U第一行與A第一行相同。4. LU分解的具體過程和順序如下：(1)第一步分解：U11 = a11(2)第二步分解：121 = a?1 / u11u12 二 a12U22 = a22 -l

5、 21U12(3)第三步分解：l 31 = a31 / U11丨32 =(a32 -l 31U12 ) / U22U13 二 a13U23 =a23 -l21U13U33 =a33 -l31 U13 - l32U23（n）第 n 步分解：依次計(jì)算：ln1、ln2Inn 4, Uin unn11_unU12Um 1I211XU22U2n1Jn1l n2.11unn 一1.3算法程序fun ctio n L,U,x=Lu_x(A,d)n, m=size(A);if n=merror(The rows and colu mns of matrix A must be equal!); return;

6、endfor ii=1: nfor i=1:iifor j=1:iiAA(i,j)=A(i,j);endendif (det(AA)=0)error(The matrix can not be divided by LU!) return;endendAn,n =size(A);L=zeros( n,n);U=zeros( n,n);for i=1: nL(i,i)=1;endfor k=1: nfor j=k: nU(k,j)=A(k,j)-sum(L (k,1:k-1).*U(1:k-1,j);endfor i=k+1: nL(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(

7、1:k-1,k)/U(k,k);endendy(1)=d(1);for i=2: nfor j=1:i-1d(i)=d(i)-L(i,j)*y(j);endy(i)=d(i)； end x(n )=y( n)/U( n,n);for i=( n-1):-1:1 for j=n :-1:i+1 y(i)=y(i)-U(i,j)*x(j);endx(i)=y(i)/U(i,i);end14實(shí)例及結(jié)果分析3xi +5x2 +6x3 +2x4 =189xi +8x2 +7x3 +5x4 =102x17x2 x3 6x4 =163x1 12x2 5x3 9x4 =17MATLAB命令窗體輸入如下： A=

8、3 5 6 2;9 8 7 5;2 7 1 6;3 12 5 9; b=18 10 16 17; L,U,x=Lu_x(A,b)得到的結(jié)果如下：A =356298752 7163 1259L =1.0000 0 03.00001.000000.6667-0.52381.00001.0000-1.00001.36960001.00003.00005.00006.00000-7.0000-11.000000-8.76190 0 02.0000-1.00004.14290.32614.7500 54.2500 -25.2500 -58.0000實(shí)驗(yàn)二Lagrange插值法數(shù)值求解2.1實(shí)驗(yàn)?zāi)康呐c要求

9、1. 了解拉格朗日插值的基本概念；2. 了解插值公式的基本原理，運(yùn)用插值公式求解實(shí)際問題；3. 編寫matlab程序，實(shí)現(xiàn)拉格朗日插值法，觀察所得結(jié)果的精確性。2.2實(shí)驗(yàn)基本原理設(shè)已知X。， xi， X2，,Xn及yi=f( xj(i=0,1,.,n), Ln(x)為不超過n次多項(xiàng)式且滿足 Ln(Xi)二 (i=0,1,.n).易知 Ln(x) =lo(x)y. ln(x)yn其中，li(x)均為n次多項(xiàng)式，再由xj (j = i)為n次多項(xiàng)式h(x)的n個(gè)根知nli(x) =C| 丨 X -Xj .最后，由j =9nli(Xj)二 C| 丨(Xi-Xj) =1 =j=0j-i,i=0,1,.

10、, n.丨丨(Xi - Xj)j =0總之,nLn(X)=二.li(x)yi，i =0nli(x)=i【j=0j朮X -XjXi _Xj.式為n階Lagrange插值公式，其中，li(x) (i=0,1,.n)稱為n階Lagrange插值的基函數(shù)。2.3算法程序fun ctio n yy=n alagr(x,y,xx)m=le ngth(x);n=len gth(y);if m=n,error(向量x與y的長度必須一致);end s=0;for i=1: nt=on es(1,le ngth(xx);for j=1: nif j=it=t.*(xx-x(j)/(x(i)-x(j);endend

11、 s=s+t*y(i);endyy=s；2.4實(shí)例及結(jié)果分析給出f(x) =1 nx的數(shù)值表，如下圖所示，用Lagrange插值計(jì)算ln(0.54)的近似值x0.40.50.60.70.8ln (x)-0.91629-0.69315-0.51083-0.35777-0.22314MATLAB命令窗體輸入如下： x=0.4 0.5 0.6 0.7 0.8; y=-0.91629 -0.69315 -0.51083 -0.35777 -0.22314 xx=0.54; n alagr(x,y,xx)得到的結(jié)果如下：ans =-0.6160用Lagrange插值法求得ln(0.54)的近似值為-0.

12、6160,與計(jì)算機(jī)求得的值-0.6162,比較，是極為相近的。實(shí)驗(yàn)三龍貝格求積公式求數(shù)值積分31實(shí)驗(yàn)?zāi)康呐c要求1. 掌握龍貝格求積的基本思路和步驟;2. 培養(yǎng)編程與上機(jī)調(diào)試能力。3.2實(shí)驗(yàn)基本原理(1)取k =0，h =b -a,求T0(0)多f(a) f(b)令1“ k ( k記為區(qū)間a,b 1的求梯形值T(寧)，即按遞推公式二分次數(shù))。h n Jf(x J 計(jì)算 T0(k)。2心k 24m 1(3) 求加速值，按公式 時(shí))= 樣羅-一m樣：(k=1,2)逐個(gè)求4 -14 -1出T表的第k行其余各元素T/kJ) ( j=1,2)。(4) 若丁 -Tk號(hào)勺(預(yù)先給出的精度)，則終止計(jì)算，并取T

13、k(0) J ；否則 令k 1 k轉(zhuǎn)(2繼續(xù)計(jì)算)。T0(0)TT1(0)T0(2)T1t2(0)TT1(2)T2T3W - -(i )Tk排成三角數(shù)表3.3算法程序fun cti on R=romberg(f, a, b, e)%參數(shù)介紹：%f -被積函數(shù)f(x)%a - x的左區(qū)間.%b -x的右區(qū)間%e -誤差限.%結(jié)果：%R-返回Romberg表n = 1; %區(qū)間二分次數(shù)while 1 %在此僅代表多次二分，在后面判斷循環(huán)終止R = zeros(n + 1, n + 1);%生成(n+1)*(n+1)的 0 矩陣R(0+1,0+1) = (b - a) / 2 * (feval(f,

14、 a) + feval(f, b); % 初始值(2 點(diǎn)梯形公式). for i = 1 : n %按照公式計(jì)算Romberg表的第一列.h = (b - a) / 2Ai;s = 0;for k = 1 : 2A(i-1)s = s + feval(f, a + (2*k - 1)*h);endR(i+1,0+1) = R(i-1+1,0+1)/2 + h*s;endfor j = 1 : n % 計(jì)算Romberg表的其他列.fac = 1 / (4Aj - 1);for m = j : nR(m+1, j+1) = R(m+1, j-1+1) + fac*(R(m+1, j-1+1) -

15、 R(m-1+1, j-1+1); endendif abs(R( n,n) - R( n+1, n+1) fun=i nli ne(4丿(l+xT),; romberg(fu n, 0,1,1e-6)得到的結(jié)果如下：t =3.1416ans =3.000003.13123.14163.14210003.13903.14163.14163.1416003.14093.14163.14163.14163.141603.10003.1333003.14143.14163.14163.14163.14163.1416實(shí)驗(yàn)四用Runge-Kutta方法求常微分方程數(shù)值解4.1實(shí)驗(yàn)?zāi)康呐c要求(1) 掌握

16、Runge-Kutta方法的基本原理；(2) 用 matlab 編寫 4 階 Runge-Kutta 程序；(3) 用寫好的程序求解常微分方程，并與系統(tǒng)自帶的ode45函數(shù)進(jìn)行比較4.2實(shí)驗(yàn)基本原理對(duì)于一階精度的歐拉公式有:yi 1 r h K1Ki 二 f(Xi,yJ當(dāng)用點(diǎn)Xi處的斜率近似值Ki與右端點(diǎn)Xi i處的斜率K2的算術(shù)平均值作為平 均斜率K*的近似值，那么就會(huì)得到二階精度的改進(jìn)歐拉公式：% i 二 yi h (Ki 心)/2Ki 二 f(Xi,yi)K2 = f (Xih,yi h Ki)依次類推，如果在區(qū)間Xi，Xii 1內(nèi)多預(yù)估幾個(gè)點(diǎn)上的斜率值Ki、K2.Km , 并用他們的加

17、權(quán)平均數(shù)作為平均斜率K*的近似值，顯然能構(gòu)造出具有很高精度的高階計(jì)算公式。經(jīng)數(shù)學(xué)推導(dǎo)、求解，可以得出四階龍格-庫塔公式，也就是在 工程中應(yīng)用廣泛的經(jīng)典龍格-庫塔算法：yn卄yn+(Ki+2K2+2K3+K4)Ki = f (Xn , yn)h hK2 = f(Xn +二，yn +：Ki)22K3 = f (Xn +, yn + K2)22I K4 = f(Xn+h,yn+hK3)4.3算法程序建立Igkt4j.m文件；fun ctio n y,x=lgkt4j(x0,x n,yO,h)x=xO:h:x n;n=len gth(x);yi=x;yi(i)=y0; for i=i: n-iKi=

18、f(x(i),yi(i);K2=f(x(i)+h/2,yi(i)+h/2*Ki);K3=f(x(i)+h/2,yi(i)+h/2*K2);K4=f(x(i)+h,yi(i)+h*K3);yi(i+i)=yi(i)+h/6*(Ki+2*K2+2*K3+K4);endy=yi;4.4實(shí)例及結(jié)果分析求微分方程 史-2xy2， y(0) = i在0 ： x : i.2時(shí)的數(shù)值解。dx(1) 首先建立f.m文件；fun ctio n z=f(x,y) z=-2*x*yA2;(2) 在MATLAB命令窗體輸入如下； x0=0;x n=1.2;y0=1;h=0.1; y,x=lgkt4j(x0,x n,y0

19、,h);n=len gth(x);fprintf(自編的龍格庫塔函數(shù)效果n); fprin tf( ix(i)y(i)n);for i=1: nfprintf(%2d %12.4f %12.4fn,i,x(i),y(i); end plot(x,y,ro-); x1,y1=ode45(f,0,1,1);fprintf(ode45 效果 n);fprin tf( ix1(i)y1(i)n);for i=1: nfprin tf(%2d %12.4f %12.4fn,i,x(i),y(i); end plot(x1,y1,bv-);運(yùn)行結(jié)果為：自編的龍格庫塔函數(shù)效果ix(i)y(i)10.0000

20、1.000020.10000.990130.20000.961540.30000.917450.40000.862160.50000.800070.60000.735380.70000.671190.80000.6098100.90000.5525111.00000.5000121.10000.4525131.20000.4098ode45效果ix1(i)y1(i)10.00001.000020.10000.990130.20000.961540.30000.917450.40000.862160.50000.800070.60000.735380.70000.671190.80000.609

21、8100.90000.5525111.00000.5000121.10000.4525131.20000.4098上文分別用自編的龍格貝塔函數(shù)和 matlab自帶的ode45函數(shù)求常微分函數(shù) 的數(shù)值解，從數(shù)值解和圖對(duì)比看出，用自帶的龍格庫塔函數(shù)求的值是非常接近的。實(shí)驗(yàn)五線性規(guī)劃任務(wù)分配問題5.1問題提出某車間有甲、乙兩臺(tái)機(jī)床，可用于加工三種工件假定這兩臺(tái)車床的可用臺(tái)時(shí)數(shù)分別為800和900,三種工件的數(shù)量分別為400、600和500，且已知用三種 不同車床加工單位數(shù)量不同工件所需的臺(tái)時(shí)數(shù)和加工費(fèi)用如下表.問怎樣分配車床的加工任務(wù)，才能既滿足加工工件的要求，又使加工費(fèi)用最低？車床 類型單位工件所

22、需加工臺(tái)時(shí)數(shù)單位工件的加工費(fèi)用可用臺(tái)時(shí)數(shù)工件1工件2工件3工件1工件2工件3甲0.41.11.013910800乙0.51.21.3111289005.2模型分析設(shè)在甲車床上加工工件1、2、3的數(shù)量分別為為、X2、X3，在乙車床上加 工工件1、2、3的數(shù)量分別為Xt、x5、x6約束條件如下:(1) 3種工件的數(shù)量分別為400、600、500,則有約束條件Xi X2 = 400、x2 x5 =600、x3 x6 =500 ；(2) 甲車床的使用時(shí)間必須小于 800,則0.4x1 1.1x2 1.0x3乞800 ；(3) 乙車床的使用時(shí)間必須小于 900,則0.5x4 1.2x5 1.3冷_900

23、 ；(4) 非負(fù)約束:務(wù) _0,i =1,2,6。此題為線性規(guī)劃模型，可用到 matlab中的優(yōu)化工具箱解線性規(guī)劃； 模型：min z = cXs.t AX MbAeq X = beq命令：x=linprog （c，A， b，Aeq,beq, VLB，VUB）5.3模型建立min 13x! 9x2 10x3 11x4 12x5 8xgX +x4 =400x2 + 花=600j X3 + & = 500 s t斗0.4/ +1.1x2 +x3 蘭 8000.5x4+1.2x5 +1.3x6 蘭900Xj 一 0,i =1,2,65.4模型求解1 0 0 1 0 0、1400 0 10 0 10X =6002 0 10 0 100丿X2 IX3x4X5改寫為：minz = 13 9 10 11 12 8 X0.4 1.1 1000 :乞00、s.tX f = 13 9 10 11 12 8; A =0.4 1.1 1 0 0 0;0 0 0 0.5 1.2 1.3;b = 800; 900;Aeq=1 0 0 1 0 0;0 1 0 0 1 0;0 0 1 0 0 1;beq=400 600 500; vlb = zeros(6,1); vub=; x,fval