《應用計算方法教程》大作業(yè).doc

應用計算方法教程作業(yè) 姓名： * * 學號： S2017* 專業(yè)： * * * * 學院： 機械工程學院 聯(lián)系方式： * 任課教師： 丁軍 作業(yè)一.用任意方法求解七階方程-5+22x-12x2+x3+13x4+10x5-7x6-4x7牛頓迭代法:y=(x)4*x7+7*x6-10*x5-13*x4-x3+12*x2-22*x+5;s=(m)28*m6+42*m5-50*m4-52*m3-3*m2+24*m-22;p0=0;N=1000;Tol=0.001;for k=1:N p1=p0-y(p0)/s(p0); if abs(p1-p0)Tol p0=p1; else break endenddisp(p1);disp(k)作業(yè)二.實驗4-2 合理利用LU分解法（使用列選主元可恰當加分）或追趕法求解下列方程組：1.3x1-x2-x3+2.4x4-3.4x5=4.22.4x1-x2-x3+1.4x4-3.7x5=6.32x1-x2-2x3+3.6x4+6.8x5=5.52.5x1-x2+4x3+3x4+6.6x5=3.61.5x1-x2-x3+5.3x4+2.8x5=6.2程序：A=1.3,-1,-1,2.4,-3.4;2.4,-1,-1,1.4,3.7;2,1,-2,3.6,6.8;2.5,-1,4,3,6.6;1.5,-1,-1,5.3,2.8;b=4.2;6.3;5.5;3.6;6.2;L=eye(5);U(1,:)=A(1,:);y=0;0;0;x=0;0;0;for k=2:5 if U(k-1,k-1)=0 disp(分解失敗);Return end L(k:5,k-1)=A(k:5,k-1)/U(k-1,k-1) U(k,k:5)=A(k,k:5)-L(k,1:k-1)*U(1:k-1,k:5) if k5 A(k+1:5,k)=A(k+1:5,k)-L(k+1:5,1:k-1)*U(1:k-1,k); endendL1=inv(L);y=L1*b;U1=inv(U);x=U1*y;disp(L);disp(U);disp(x);結(jié)果：作業(yè)三.實驗4-311/21/31/41/51/61/71/81/91/101/21/31/41/51/61/71/81/91/101/111/31/41/51/61/71/81/91/101/111/121/41/51/61/71/81/91/101/111/121/131/51/61/71/81/91/101/111/121/131/141/61/71/81/91/101/111/121/131/141/151/71/81/91/101/111/121/131/141/151/161/81/91/101/111/121/131/141/151/161/171/91/101/111/121/131/141/151/161/171/181/101/111/121/131/141/151/161/171/181/19觀察結(jié)論：矩陣的條件數(shù)越來越大，而且遠遠大于1，這是一個病態(tài)問題 1） LU分解法算法實現(xiàn)：% n=5:10clear;clc;n=10for i=1:n for j=1:n H(i,j)=1/(i+j-1); endenddisp(H)cond(H,1)format short;n=size(H);L=eye(n);U(1,:) = H(1,:); for k = 2:n if U(k-1,k-1)=0 disp(分解失敗);return end L(k:n,k-1)=H(k:n,k-1)/U(k-1,k-1); U(k,k:n)=H(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n); if kn H(k+1:n,k)=H(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k); end end B(1:10,1)=1; Y=LB; L U X=UY 運行結(jié)果：L=1.00000000000000.50001.0000000000000.33331.00001.000000000000.25000.90001.50001.00000000000.20000.80001.71432.00001.0000000000.16670.71431.78572.77782.50001.000000000.14290.64291.78573.33334.09093.00001.00000000.12500.58331.75003.71215.56825.65383.50001.0000000.11110.53331.69703.95966.85318.61547.46674.00001.000000.10000.49091.63644.11197.930111.630812.60009.52944.50001.0000U=1.00000.50000.33330.25000.20000.16670.14290.12500.11110.100000.08330.08330.07500.06670.05950.05360.04860.04440.0409000.00560.00830.00950.00990.00990.00970.00940.00910000.00040.00070.00100.00120.00130.00140.001500000.00000.00010.00010.00010.00020.0002000000.00000.00000.00000.00000.00000000000.00000.00000.00000.000000000000.00000.00000.0000000000000.00000.00000000000000.0000X=1.0e+06 *-0.00000.0010-0.02380.2402-1.26113.7834-6.72617.0007-3.93790.92372） LU分解迭代求精法算法實現(xiàn)：% n=5:10clear;clc;n=10;for i=1:n for j=1:n H(i,j)=1/(i+j-1); endendformat short;n=size(H);L=eye(n);U(1,:) = H(1,:);tol=10(-3); for k = 2:n if U(k-1,k-1)=0 disp(分解失敗);return end L(k:n,k-1)=H(k:n,k-1)/U(k-1,k-1); U(k,k:n)=H(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n); if kn H(k+1:n,k)=H(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k); end end B(1:10,1)=1; Y=LB; X=UY; r=B-H*X; Y=Lr; w=UY; format long while norm(w)tol r=B-H*X; Y=Lr; w=UY; X=X+w; end X運算結(jié)果：X =1.0e+06 *-0.0000099983463550.000989857965323-0.0237569902910640.240212759410420-1.2611305603384353.783425339621443-6.7261397788656957.000720697677755-3.9379270271364350.923715699219681 結(jié)果分析：對于同一復雜程度的方程組，LU分解迭代求精法所得的解較LU分解法所得到的解的精確度更高。作業(yè)四.實驗5-1 分別用Jacobi、Seidel、Sor（=0.8,1.1,1.2,1.3,1.4,1.5）迭代求解下面的方程組，并做結(jié)果分析。初值x(0)=(0,0,0,0,0)T，精度要求：x(k+1)-xkx(k+1)10-5。12.3x1-2x2-x3+3.4x4-3.7x5=41.4x1+9x2-3x3+2.4x4+2.7x5=2.32.1x1+x2+8x3+2.6x4+5.8x5=2.53.5x1-2.1x2+x3+13x4+4.6x5=3.62.5x1-x2-2x3+5.3x4+14.8x5=2.2Jacobi:A=12.3,-2,-1,3.4,-3.7;1.4,9,-3,2.4,2.7;2.1,1,8,2.6,5.8;3.5,-2.1,1,13,4.6;2.5,-1,-2,5.3,14.8;b=4.8;2.3;2.5;3.6;2.2;X0=0;0;0;0;0;X=X0;K=1;while K=30for i=1:5X(i)=(b(i)-A(i,:)*X0)/A(i,i)+X0(i);endif norm(X-X0)/X0.00001disp(X);disp(K);return;endK=K+1;X0=X;end結(jié)果：Seidel:A=12.3,-2,-1,3.4,-3.7;1.4,9,-3,2.4,2.7;2.1,1,8,2.6,5.8;3.5,-2.1,1,13,4.6;2.5,-1,-2,5.3,14.8;b=4.8;2.3;2.5;3.6;2.2;X0=0;0;0;0;0;X=X0;K=1;while K=20 for i=1:5 X(i)=(b(i)-A(i,:)*X)/A(i,i)+X(i); end if norm(X-X0)/X0.00001 disp(X); disp(K); return; end K=K+1; X0=X;enddisp(發(fā)散);結(jié)果：SOR:A=12.3,-2,-1,3.4,-3.7;1.4,9,-3,2.4,2.7;2.1,1,8,2.6,5.8;3.5,-2.1,1,13,4.6;2.5,-1,-2,5.3,14.8;b=4.8;2.3;2.5;3.6;2.2;X0=0;0;0;0;0;X=X0;K=1;w=0.8,1.1,1.2,1.3,1.4,1.5; for j=1:6 m=w(j);while K=20 for i=1:5 X(i)=m*(b(i)-A(i,:)*X)/A(i,i)+X(i); end if norm(X-X0,inf)/norm(X,inf) syms x; p0=1; p1=x-quad(x,0.40,1.05)/(1.05-0.40)*p0輸出：p1 =x - 29/40Step2：p2=x2-quad(x.2,0.40,1.05)/(1.05-0.40)*p0-quad(x.3-29*x.2/40,0.40,1.05)/quad(x-29/40).2,0.40,1.05)*p1輸出：p2 =x2 - (29*x)/20 + 1177/2400 Step3： p3=x3-quad(x.3,0.40,1.05)/(1.05-0.40)*p0-quad(x.4-29*x.3/40,0.40,1.05)/quad(x-29/40).2,0.40,1.05)*p1-quad (x.5 - 29*x.4/20 + 1177/2400*x.3,0.40,1.05)/ quad(x.2 - (29*x)/20 + 1177/2400).2,0.40,1.05)*p2輸出：p3 =x3 - (87*x2)/40 + (3027*x)/2000 - 53621/160000Step4： x=0.40,0.55,0.65,0.80,0.90,1.05;y=0.41075;0.57815;0.69675;0.88811;1.02652;1.25386;A(1:6,1)=ones(size(x);A(1:6,2)=x - 29/40 ;A(1:6,3)=x.2 - (29*x)/20 + 1177/2400 ;A(1:6,4)=x.3 - (87*x.2)/40 + (3027*x)/2000 - 53621/160000; C=Ay輸出：C = 0.8042 1.2881 0.3989 0.2137Step5： s=0;k=1;while k s輸出：s = 0.0075 結(jié)果分析：比較(1)和（2）可以發(fā)現(xiàn)正交之后得到的最小二乘法的誤差是直接多項式擬合的萬分之一，通過正交的pj(x)的函數(shù)系解決最小二乘法的病態(tài)結(jié)果問題，提供了求最小二乘法擬合多項式的快捷方法。作業(yè)七.計算 （1）編寫復化梯形公式和復化Simpson公式通用子程序，分別采用4,8,16,32,64等分區(qū)間計算。（2）使用Romberg求積公式。(1) 復化梯形公式程序：a=0;b=1;s=0;f=(x)(log(1+x)/x);fa=1;fb=log(2);n=2;for i=1:5n=n*2;h=(b-a)/(n);f1=0;for i=1:n-1f1=f1+f(a+h);ends=h*(fa+fb+2*f1)/2;disp(s)end結(jié)果：復化simposon公式程序：a=0;b=1;s=0;f=(x)(log(1+x)/x);fa=1;fb=log(2);n=1;for j=1:5 n=n*2;h=(b-a)/(2*n);f1=0;for i=0:n-1f1=f1+f(a+(2*i+1)*h);endf2=0;for i=1:n-1f2=f2+f(a+2*i*h);ends=h*(fa+fb+4*f1+2*f2)/3;disp(s)end結(jié)果：(2)程序：function y=romberg(f,n,a,b)f=(x)(log(1+x)/x);a=0;b=1;n=4;z=zeros(n,n);h=b-a;z(1,1)=(h/2)*(f(a)+f(b);f1=0;fori=2:nfork=1:2(i-2)f1=f1+f(a+(k-0.5)*h);endz(i,1)=0.5*z(i-1,1)+0.5*h*f1;h=h/2;f1=0;forj