《应用计算方法教程》大作业.doc

应用计算方法教程作业 姓名： * * 学号： S2017* 专业： * * * * 学院： 机械工程学院 联系方式： * 任课教师： 丁军 作业一.用任意方法求解七阶方程-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)作业二.实验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);结果：作业三.实验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观察结论：矩阵的条件数越来越大，而且远远大于1，这是一个病态问题 1） LU分解法算法实现：% 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 运行结果：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分解迭代求精法算法实现：% 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运算结果：X =1.0e+06 *-0.0000099983463550.000989857965323-0.0237569902910640.240212759410420-1.2611305603384353.783425339621443-6.7261397788656957.000720697677755-3.9379270271364350.923715699219681 结果分析：对于同一复杂程度的方程组，LU分解迭代求精法所得的解较LU分解法所得到的解的精确度更高。作业四.实验5-1 分别用Jacobi、Seidel、Sor（=0.8,1.1,1.2,1.3,1.4,1.5）迭代求解下面的方程组，并做结果分析。初值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结果：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(发散);结果：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 结果分析：比较(1)和（2）可以发现正交之后得到的最小二乘法的误差是直接多项式拟合的万分之一，通过正交的pj(x)的函数系解决最小二乘法的病态结果问题，提供了求最小二乘法拟合多项式的快捷方法。作业七.计算 （1）编写复化梯形公式和复化Simpson公式通用子程序，分别采用4,8,16,32,64等分区间计算。（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结果：复化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结果：(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