科学工程计算

2010届科学与工程计算作业指导教师:丁军作 者：陈永学 号：G201081372010年11 月2 日目录作业1：应用牛顿迭代法和斯蒂芬森加速法解方程 11.1牛顿迭代法： 11.2斯蒂芬森(Steffensen) 加速法2作业 2：线性方程组数值解法程序设计与验证42.1Gauss按比例歹U选主元消去法52.2Jacobi 迭代法72.3 Seidel 迭代法8作业 3： 曲线拟合93.1 牛顿插值 103.2拉格朗日插值113.3 三次拟合曲线123.4四次拟合曲线133.5四种插值拟合值图形比较14作业四：几种求积公式积分154.1 复化梯形公式求解：154.2 复化 simpson 公式求解： 164.3 龙贝格积分法计算17作业五：常微分方程初值问题的数值解法185.1改进euler方法185.2四级四阶RK方法20作业1：应用牛顿迭代法和斯蒂芬森加速法解方程【题目】使用牛顿迭代法和斯蒂芬森Steffensen)加速法求解xA3-5xA2+6x-1=0全部根，(注：根全部在区间D,4】中)要求精确到10A(-6)并给出以-5,5区间的100等分点(即-5， -4.9,-4.8,4.8,4.9,5为) 初始点，收敛到不同解的示意图1.1牛顿迭代法：算法说明：牛顿法递推公式为：f ( x )X X nn+1n f (x )(1.1)n取-5， -4.9,-4.8,4.8,4.9,5为初始点，分别代入式1.1，求出满足条件的解(由于题目中给出其根全部在区间D,4】中，此方法求出的解即为KA3-5xA2+6x-1=0全部根)源代码：tol=10A(-6)；%误差限N=2000;%最大迭代步骤i=1；for p0=-5:0.1:5;for k=1:Np1=p0-(p0A3-5*p0A2+6*p0-1)./(3*p0A2-10*p0+6);if abs(p1-p0)tolbreakendp0=p1;enda(i)=p1;i=i+1;endp0=-5:0.1:5c0=find(abs(a-0.1981)0.1),c1=find(abs(a-1.555)0.1),c2=find(abs(a-3.2470)0.001) subplot(111)plot(p0(c0),0,b-*,p0(c1),0,r-d,p0(c2),0,c-o);gridaxis(0 4 -10 10);输出结果：1.2斯蒂芬森Steffensed加速法算法说明：Steffensen加速法递推公式为：廿 + 一厂 f (x+ f(xk ) - f (x) f Wk -1k-1k -1,且有x1f (a)f (a + f (a) - f (a)f(a)源代码：tol=10人(-6);%误差限N=2000;%最大迭代步骤i=1;n=0;for p0=-5:0.1:5p(1)=p0;while n=Nfor k=1:2p(k+l)=p(k)-(p(k)人3-5*p(k)人2+6*p(k)-l)/(3*p(k)人2-10*p(k)+6牛顿 迭代法endp1=pW-(p (2)-p(1)A2/(p -2*p (2)+p(1);%Aitken 加速f abs(p1A3-5*p1A2+6*p1-1)vtolbreakendn=n+1;p(1)=p1;enda(i)=p1;i=i+1;endp0=-5:0.1:5c0=find(abs(a-0.1981)0.1),c1=find(abs(a-1.555)0.1),c2=find(abs(a-3.2470)0.001) subplot(111)plot(p0(c0),0,b-*,p0(c1),0,r-d,p0(c2),0,c-o);gridaxis(0 4 -10 10);输出结果：作业2：线性方程组数值解法程序设计与验证【实验目的及意义】1）掌握Gauss消去法及Gauss列主元消去法，能用这两种方法求解方程组。2）掌握矩阵和向量范数的定义，了解它们的性质。3）掌握Jacobi迭代法，能应用Tacobi迭代法求解方程组；4）掌握Seidel迭代法，能应厝eide迭代法求解方程组。5）能用矩阵的范数判断迭代法的收敛性。【实验要求】设计和验证Gauss按比例列选主元消去法和Jacobi迭代法、Seidel迭代法算法的程序17.0310.6152.9111.0071.0060.000 0.230 1.00034.2111.0002.1000.3001.70052.3220.0000.50013.0000.5001.0001.500V 544.5013.1103.90761.70512.1708.999x240.2360.1018.0120.0170.9104.6180.10029.304 1.000234.5521.803 ,一117.818,误差不超过 = 10-6 .2.1Gauss按比例列选主元消去法算法说明：为了克服Gauss消元法的两点不足，在每步消元时先选择绝对值较大的数作为乘数m 的分母，控制舍入误差的扩大。列主元消元法是一种局部主元消元法，按比例列选主 ik元消去法是在列主元消元法的基础上给出的一种改进，主要作用是消除方程中比例因子 对选主元的影响。源代码：A=17.031 -0.615 -2.911 1.007 -1.006 0.000;-1.000 34.211 -1.000 -2.100 0.300 - 1.700;0.000 0.500 13.000 -0.500 1.000 -1.500;4.501 3.110 -3.907 -61.705 12.170 8.999;0.101 - 8.012 -0.017 -0.910 4.618 0.100;1.000 2 3 4.5 5 21.803;b=0.230;-52.322;54;240.236;29.304;-117.818;B=A,bTol=10人(-6);n,m=size(B);if n=m-1errorreturnendfor i=1:nW(i)=max(B(i,:);endfor k=1:n-1K=B(k:n,k)./W(k:n);l,m=max(K);if m=1T=B(k,:);B(k,:)=B(k+m-1,:);B(k+m-1,:)=T;S=W(k);W(k)=W(k+m-1);W(k+m-1)=S;endfor i=k+1:1:nB(i,:)=B(i,:)-B(k,:)*B(i,k)/B(k,k);endif B(n,n)=0errorreturnelse X(n)=B(n,n+1)/B(n,n);endfor k=n-1:-1:1X(k)=(B(k,n+1)-dot(B(k,k+1:n),X(k+1:n)/B(k,k);if abs(X(k)-X(k+1)=1.0e-6x0=y; y=B*x0+f;n=n+1;endyn A=17. 031 -0. 615 -2. 911 1. 007 -1. 006 0.000;-1.000 34.211 -1.000 -2b=0. 230;-52. 322;54;240.236;29.304;-117.818;x0=zeros (1 lengt?L(b);Tol=10rt(-S);N=100;n. n=size (A);k=l ;while k=Nfor i=l:nx(i)二(b(i)-A(i, :)*K0?-Hi0(i)*A(iJ i)/A(i, i);endif abs(x-xO)Tolbre：ak;endk=k+l;x0=x;enddisp (x)0.8959-1.96403.2807-4.47932. 1625-5.28752.3 Seidel迭代法算法说明：严=2卜寻严迄帶源代码：A=17.031 -0.615 -2.911 1.007 -1.006 0.000;-1.000 34.211 -1.000 -2.100 0.300 - 1.700;0.000 0.500 13.000 -0.500 1.000 -1.500;4.501 3.110 -3.907 -61.705 12.170 8.999;0.101 - 8.012 -0.017 -0.910 4.618 0.100;1.000 2 3 4.5 5 21.803;b=0.230;-52.322;54;240.236;29.304;-117.818;x0=0;0;0;0;0;0;Tol=10人(-6);N=100;n,n=size(A);x=x0;k=1;while k=Nfor i=1:n x(i)=(b(i)-A(i,:)*x+x(i)*A(i,i)/A(i,i);endif abs(x-x0)Tolbreak;endk=k+1;x0=x;enddisp(x)输出结果： A=17. 031 -0. 615 -2. 911 1. 007 -1. 006 0. 000;-l. 000 34. 211 -1. 000 -2b= Ll. 230;-52. 322; 54; 240. 236; 29. 304;-U7. 818;x0=U;0;0;0;0;0;Tol=10A(-6);N=100;n, n=size (A);x=x0;k=l;while k=Nfor i=1:nx (i)=(b (i)-A(i,:)机知(i) *A(i, i)/A(i, i);endif abs(x-xO)Tolbre：dk;endk=k+l;x0=x;已nddisp (x)0. 8959-1. 96403. 2807-4.47932.1625-5. 2875I作业 3： 曲线拟合x=0,1,2,3,4,5,6,7,8,910y=7.579,8.454,6.819,4.620,3.161,3.108,4.491,6.698,8.482,7.956,2.594使用以下方法求0到10的100等分点(即0.1,0.2,0.3,9.8,9.9)的插值或拟合值(1) 牛顿插值(11点全用)(2) 任何点选择离其最近的四个节点做拉格朗日差值(如3点.4 使用节点2,3,4,5进行插值，8.8 使用节点7,8,9,10插值(3) 三次拟合曲线(4) 四次拟合曲线(5) 上述四种插值或拟合值在同一图上画图比较3.1牛顿插值算法说明：Nn(x) = f(x0)+ fx0Jxa(x- x0)+ fx0JxLJx2(x-x0)(x- Xj) + + fx0J -. xn(x - x0) - (x - xn_ L)源代码：x=0 1 2 3 4 5 6 7 8 9 10;y=7.579 8.454 6.819 4.620 3.161 3.108 4.491 6.698 8.482 7.956 2.594;n=10;for j=1:nfor i=n+1:-1:j+1y(i)=(y(i)-y(i-1)/(x(i)-x(i-j);endendy1=y(n+1);x1=0.1;for k=1:99for j=n:-1:1y1=y(j)+(x1-x(j)*y1;endx2(k)=x1;y2(k)=y1;x1=x1+0.1;y1=y(n+1);end3.2拉格朗日插值算法说明：首先将x=l:10的100的等分，得101个值，每一个值取整，求其最近的 4 个数来做插值，这样每一个值都得到4 个数为一组的插值，然后利于拉格朗日插值(x x )公式：Ln(x)=加n丄* y ,将每组数的4个值代入得到各自的插值基函数，并nX Xkk=0 j=0k j 丿j丰k将其对应的x1=0:0.1:10中的值代入到所求基函数中得到对应的插值，及存入到数组1中。源代码：x=0 1 2 3 4 5 6 7 8 9 10;y=7.579 8.454 6.819 4.620 3.161 3.108 4.491 6.698 8.482 7.956 2.594;x1=0:0.1:10;for k=1:101y1(k)=0;i=floor(x1(k);if i1for m=0:3t=1;for j=0:3if j=mt=t*(x1(k)-x(j+1)/(x(m+1)-x(j+1);endendy1(k)=y1(k)+t*y(m+1);endelseif ieps)J=J+1; h=h/2; x=a+h:2*h:b-h;R(J+1,1)=R(J,1)/2+h*sum(subs(f,x,x);for k=1:min(3,J) R(J+l,k+l)=R(J+l,k)+(R(J+l,k)-R(J,k)/(4人k-1);end if(J3)err=abs(R(J+l,4)-R(J,4);endendquad=R(J,4)输出结果： syms x;a=0;b=4;h=b-a;eps=10(-5);err=l;J=0;R=zeros(4, 4);f=exp(x);f a=subs (fj ?a);fb=subs(fj ? j b);Rd, l)=h*(fa+ft.)/2;whi1e(erreps)J=J+1； h=h/2;x=a+h:2*h:b-h;R(J+lj l)=R(Jj l)/2+h*sujrL(ubs (f/i-f x);for k=l: min(3j J)RQ+lJk+l)=R(J+lJk) + (R(J+lJk)-R(JJk)/(4Ak-l);endif (J3)err=abs (R (J+lj 4) -R (L 4);endendquad=R (Jj 4)quad =作业五：常微分方程初值问题的数值解法分别使用改进euler方法和四级四阶RK方法求解10-3题取步长h=0.01,计算到y(0.5)5.1改进euler方法算法说明：%+：=% + hfOtpyJ如+丄=% +f%+ f(Xi+i屁+源代码：h=0.01;n=50;x=0;y=0;for i=1:nfl=-x人2+x-y;y1=y+h*f1;f2=-(x+h)人2+x+h-yl;y2=y+h*f2;y=(yl+y2)/2;a(i)=y;x=x+h;enddisp(a);输出结果：Columns 1 through 60.00000.00020.00040.00080.00120.0017Columns7 through :120.00230.00290.00370.00450. 00540.0064Columns13 through180.00740.00850.00960.01080.01210. 0134Columns19 through240.01480.01620.01770.01920.02070.0223Columns25 through300.02390.02560.02720.02890.03070.0325Columns 31 through 360.03420.03600.03790.03970. 04160. 0434Columns 37 through 420.04530. 04720.04910.05100.05280. 0547Columns 43 through 480.05660.05850.06040.06220. 06410.0659Coluntns 49 through 500.06780.06965.2四级四阶RK方法算法说明：对于微分方程，经典龙格库塔法的计算公式为K = f (x , y )1K2 K3K4yn+1nhh=f (x + 2, y + 2 K1)n 2 n 2 1hh=f (x + ，y + 2 K 2)2 2 2=f (x + h, y + hK )3h=y + (K + 2 K + 2 K + K )6134源代码：x=0;y=0;n=50;h=0.01;for i=1:nf1=-x人2+x-y;k1=h*f1;f2=-(x+h/2)人2+x+h/2-y-k1/2; k2=h*f2;f3=-(x+h/2)A2+x+h/2-y-k2/2; k3=h*f3;f4=-(x+h)A2+x+h-y-k3; k4=h*f4;y=y+(k1+2*k2+2*k3+k4)/6;b(i)=y; x=x+h;end disp(b); 输出结果：0. i：ii：ii：ii：i0. I：II：II：I20. I：IOI：I40. COOS0.00120.00170.0023ColujTLiLS 15 through280.00960.01080. 01210. 01340. 01480. 01620. 0177CoIi-lultls 29 through 420.03070.03250. 03420. 03600. 03790. 03970. 0416ColujTLiLS 43 through500.05660.05850. 06040. 06230. 06410. 06600.0678ColujTLns 1 thr cmgh 140. 06960. 00290. 00370.00450.00540.00640.00740.00850. 01920.02070.02230. 02390.02560.02720.02900. 04340. 04530. 04720. 04910.05100.05290.0547