数值分析作业.doc

实验报告一题目：分析误差的方法与原则摘要：用计算机解决实际问题就要先建立数学模型，而模型与实际问题之间会出现误差无可避免。因此在数值分析中除了要研究数学问题的数值方法与理论，还要研究计算结果的误差是否满足精度要求，这就是误差估计问题。 前言：（目的和意义）掌握截断误差与舍入误差。数学原理：截断误差指数学模型与数值方法之间的误差.由于实际运算只能完成有限步，因此就要用极限的方法对无限过程进行简化，这样产生的误差成为截断误差。舍入误差指的是由于计算机的位数有限带来的误差称为“舍入误差”，计算机数字位数有限，对数据往往要进行四舍五入取其近似值从而产生误差。舍入误差是为了计算的方便人为造成的，而截断误差是由于计算方法造成的不可避免的误差。1.绝对误差：，简称误差。2. 相对误差：，在实际的计算过程中间，通常取3. 相对误差限： （）4. （1）计算In=exp(-1)x n exp(x)dx x在（0,1）中 (n=0,1,)并估计误差。解：由分部积分得In的递推公式 In=1-nIn-1（n=1,2，）,程序如下：I=int(x0*exp(x-1),x,0,1)I = 1 - 1/exp(1) vpa(I,4) ans = 0.6321 I=int(x1*exp(x-1),x,0,1) I = 1/exp(1) vpa(I,4) ans = 0.3679 I=int(x2*exp(x-1),x,0,1) I = 1 - 2/exp(1) vpa(I,4) ans = 0.2642 I=int(x3*exp(x-1),x,0,1) I = 6/exp(1) - 2 vpa(I,4) ans = 0.2073I=int(x4*exp(x-1),x,0,1) I = 9 - 24/exp(1) vpa(I,4) ans = 0.1709 I=int(x5*exp(x-1),x,0,1) I = 120/exp(1) - 44 vpa(I,4) ans = 0.1455 I=int(x6*exp(x-1),x,0,1) I = 265 - 720/exp(1) vpa(I,4) ans = 0.1268 I=int(x7*exp(x-1),x,0,1) I = 5040/exp(1) - 1854 vpa(I,4) ans = 0.1124 I=int(x8*exp(x-1),x,0,1) I = 14833 - 40320/exp(1) vpa(I,4) ans = 0.1009 I=int(x9*exp(x-1),x,0,1) I = 362880/exp(1) - 133496 vpa(I,4) ans = 0.09161此题迭代到第八步。例14要使20的近似值的相对误差限小于0.1%，要取几位有效数字？解：近似数x*，若x*具有n位有效数字，则其相对误差限为，*(1/2a1)*10-(n-1) 则其至少有n位有效数字。此题中20的近似值=4.4知a1=4，故要取n=4.就有20的近似值的相对误差限小于0.1%。 其编程如下：I=sym(20);sqrt(I)ans = 2*5(1/2) vpa(sqrt(I),4) ans = 4.472（3）a0=3,ak=2ak-1,用秦九韶算法求p100(0.5),p150(13)掌握秦九韶算法的基本原理与应用。5.秦九韶算法：.程序设计：秦九韶算法的主程序：function Q=Qjs(n,x)a(1)=3;for i=1:n+1 a(i+1)=2.*a(i)+3;endS1=a(n+1);for j=1:n S=x*S1+a(n+1-j); S1=S;endQ=S;end结果分析和讨论：1 已知Pnx=anxn+an-1x(n-1)+.+a1x+a0, an=2 an-1+3,a0=3,求P100(0.5),P150(13).解：计算结果为Qjs(100,0.5)ans = 600Qjs(150,13)ans = 1.0995e+213结论：多项式的求解有多种计算方法，秦九韶算法是一种将一元n次多项式的求值问题转化为n个一次式的算法，一般，一元n次多项式的求值需要经过（n*（n+1）/2次乘法和n次加法，而秦九韶算法只需要n次乘法和n次加法。 实验报告二题目： 拉格朗日插值法和牛顿插值法求近视值摘要：许多实际问题都用函数y=f（x）来表示某种内在规律的数量关系，其中一部分函数是通过实验或观测得到的。很多f（x）只能给出a,b上的一系列点xi的函数值或者一张函数表，有的函数虽然有解析表达式，但由于计算复杂，我们希望可以根据给定的函数表做一个既能反映函数f(x)的特性，又便于计算的简单函数P(x)，用P(x)近似f(X)。通常选一类简单的函数作为P(x)，并使P(xi)=f(xi)对i=1,2,n成立。由此确定的P(x)就是我们希望的插值函数，此即为插值法。前言：（目的和意义）掌握Lagrange插值算法。数学原理：拉格朗日插值法n次拉格朗日插值多项式为：Ln(x)=y0l0(x)+y1l1(x)+y2l2(x)+ynln(x)n=1时，称为线性插值，L1(x)=y0(x-x1)/(x0-x1)+ y1(x-x0)/(x1-x0)=y0+(y1-x0)(x-x0)/(x1-x0)n=2时，称为二次插值或抛物线插值，精度相对较高L2(x)=y0(x-x1)(x-x2)/(x0-x1)/(x0-x2)+y1(x-x0)(x-x2)/(x1-x0)/(x1-x2)+y2(x-x0)(x-x1)/(x2-x0)/(x2-x1) 牛顿插值法：现构造一个插值多项式Nn(x)，利用牛顿插值公式，增加一个节点时，只需在后面多计算一项，而前面的计算仍有用；另一方面Nn(x)的各项系数恰好又是各阶差商，而各阶差商可用差商公式来计算。由线性代数知，对任何一个不高n次的多项式P(x)=b0b1xb2x2bnxn (幂基) 也可将其写成P(x)=a0a1(xx0)a2(xx0) (xx1)an(xx0) (xxn-1)其中ai为系数，xi为给定节点，可由求出ai，一般情况下，牛顿插值多项式Nn(x)可写成：Nn(x)= a0a1(xx0)a2(xx0) (xx1)an(xx0) (xxn-1)线性插值法： 线性插值法：线性插值经常用于已知函数 f 在两点的值要近似获得其它点数值的方法，这种近似方法的误差定义为： R（T）=f（x）-P（x）其中 p 表示上面定义的线性插值多项式 根据罗尔定理，我们可以证明：如果 f 有二阶连续导数，那么误差范围是 函数上两点之间的近似随着近似的函数的二阶导数的增大逐渐变差。直观上来看：函数的曲率越大，线性插值近似的误差也越大。 二、主要思路使用线性方程组求系数构造插值公式相对复杂，可改用构造方法来插值。对节点xi(i=0,1,n)中任一点xk(0=k=n)作一n 次多项式lk(xk)，使它在该点上取值为1，而在其余点xi(i=0,1,k-1,k+1,n)上为0，则插值多项式为Ln(x)=y0l0(x)+y1l1(x)+y2l2(x)+ynln(x)上式表明：n 个点xi(i=0,1,k-1,k+1,n)都是lk(x)的零点。求得lk三计算方法及过程：（拉格朗日插值法）1.输入节点的个数n 2.输入各个节点的横纵坐标 3.输入插值点 4.调用函数，返回（牛顿插值法）1.先后输入节点个数和节点的横纵坐标，插值点的横坐标，最后输入精度。2. 用do-while循环语句得到跳出循环时k的值3.将k值与n-1进行比较，若达到精度时k x0=0.4 0.5 0.6 0.7 0.8; y0=-0.916291 -0.693147 -0.510826 -0.356675 -0.223144;format long interp1(x0,y0,0.54)ans = -0.620218600000000一次 x0=0.4 0.5 0.6 0.7 0.8;y0=-0.916291 -0.693147 -0.510826 -0.356675 -0.223144;x=0.54; n=length(x0);s=0; for j=0:(n-1) t=1; for i=0:(n-1) if i=j t=t*(x-x0(i+1)/(x0(j+1)-x0(i+1); end end s=s+t*y0(j+1);ends s = -0.616142715200000二次 interp1(x0,y0,0.54,spline)ans = -0.6160998500000021：解：程序如下：x0=-5:5;y0=1./(1+x0.2);x=-5:0.05:5;m=length(x);for k=1:mtx=x(k); n=length(x0); s=0; for j=0:(n-1) t=1; for i=0:(n-1) if i=j t=t*(tx-x0(i+1)/(x0(j+1)-x0(i+1); end end s=s+t*y0(j+1); end yx(k)=s;endans=Ih(x)=(x(i+1)-x)*(1/(1+x(i)*x(i)+(x-x(i)*(1/(1+1+x(i+1)*x(i+1)ans=x=-4.5|4.5 f(x)=0.0471 Ih(x)=0.0471x=-3.5|3.5 f(x)=0.0755 Ih(x)=0.0794x=-2.5|2.5 f(x)=0.1379 Ih(x)=0.1500x=-1.5|1.5 f(x)=0.3077 Ih(x)=0.3500x=-0.5|0.5 f(x)=0.800 Ih(x0=0.7500ans=max abs(f(x0-Ih(x) for i=1:n-1h(i)=x(i+1)-(i);endfor i=1:n-2 k(i)=h(i+1)/(h(i)+h(i+1); u(i)=h(i)/(h(i)+h(i+1);endfor i=1:n-2 al(i)=3*(u(i)*(y(i+2)-y(i+1)/h(i+1)+k(i)*(y(i+1)-y(i)/h(i);endg0=3*(y(2)-y(1)/h(n-1);g00=3*(y(n)-y(n-1)/h(n-1);g=g0,g1,g00;g=ranspose(g);k1=k 1;ui=i u;Q=2*eye(s)+diag(u1,1)=diag(x1,-1);m=ranspose(Q/q);syms x;p1(i)=(1+2*(x-x(i)/h(i)*(x-x(i+1)/h(i)2*y(i);p2(i)=(1-2*x(i+1)/h(i)*(x-x(i)/h(i)2*y(i+1);p3(i)=(x-x(i)*(x-x(i+1)/h(i)2*m(I);p4(i)=(x-x(i+1)*(x-x(i)/h(i)2*m(i);p(i)=p1(i0+p2(i)+p3(i)+p4(i);p(i)=simple(p(i);ends1=p(1);s2=p(2);s3=p(3);s4=p(4)ans=s(1)=-6.2697*(x-0,25)3+10*(0.3-x)+10.9697*(x-0.25)s(2)=-3.4831*(o,39-x)3-1.5956*(x-0.3)3+6.1138*(0.39-x)s(3)=-2.3933*(0.45-x)3-2.8622*(x-0.39)2+10.4169*(0.45-x)+11.1903*(x-0.39)s(4)=-2.1467*(0.53-x)3+8.3987*(0.53-x)+9.1*(x-0.45)结论：各种插值法都有自己利与弊，拉格朗日插值法运算相对复杂，但组成抛物差值，精度可以提高很多.牛顿插值法及拉格朗日插值法等线性插值法只适合在已知点较少的情况下使用，当坐标点很多时应将区间分成小段进行分段差值或是分段抛物差值。掌握多个思想方法组合的技巧 实验报告三一：题目：最小二乘曲线拟合二：前言：（目的和意义）通过上机计算，对曲线的最小二乘法的拟合有进一步的掌握，并且能够熟练的运用这种方法三，实验原理：在实验数据处理中，往往要根据一组给定的实验数据，求出自变量x与因变量y的函数关系，这是为待定参数，由于观测数据总有误差，且待定参数ai的数量比给定数据点的数量少(即nm)，因此它不同于插值问题.这类问题不要求通过点，而只要求在给定点上的误差的平方和最小.当时，即(4.4.1)这里是线性无关的函数族，假定在上给出一组数据，以及对应的一组权，这里为权系数，要求使最小，其中这就是最小二乘逼近，得到拟合曲线为y=s(x)，这种方法称为曲线拟合的最小二乘法.(4.4.2)中实际上是关于的多元函数，求I的最小值就是求多元函数I的极值，由极值必要条件，可得(4.4.3)根据内积定引入相应带权内积记号(4.4.4)则(4.4.3)可改写为这是关于参数的线性方程组，用矩阵表示为(4.4.5)(4.4.5)称为法方程.当线性无关，且在点集上至多只有n个不同零点，则称在X上满足Haar条件，此时(4.4.5)的解存在唯一(证明见3).记(4.4.5)的解为 从而得到最小二乘拟合曲线(4.4.6)可以证明对，有故(4.4.6)得到的即为所求的最小二乘解.它的平方误差为(4.4.7)均方误差为在最小二乘逼近中，若取，则，表示为最小二乘的二次拟合效果如下:23观察物体的直线运动，得出以下数据 ，求运动方程时间t(秒)00.91.93.03.9 5.0距离（米）010305080 110求其拟合多项式.解：function c=mypolyfit(x,y,m)n=length(x);b=zeros(1,m+1);f=zeros(n,m+1);for k=1:m+1 f(:,k)=x.(k-1);enda=f*f;b=f*y;c=ab;c=flipud(c);另存为mypolyfit.m后,以下进入拟合阶x=0 0.9 1.9 3.0 3.9 5.0;y=0 10 30 50 80 110;m=1;C=mypolyfit(x,y,m)C =22.2538 -7.8550 plot(x,y,ok,LineWidth,2);grid on;title(最小二乘一次拟合效果图);x0=0:0.5:5;y0=C(1).*x0+C(2);hold on ,plot(x0,y0,-r);最小二乘的一次拟合效果如下: m=2; C=mypolyfit(x,y,m)C = 2.2488 11.0814 -0.5834 figure; plot(x,y,ok,LineWidth,2); grid on; title(最小二乘二次拟合效果图); x0=0:0.5:5; y0=C(1).*x0.2+C(2).*x0+C(3); hold on ,plot(x0,y0,-r); 24在某化学反应里，根据实验所得分解物浓度与时间关系如下：时间t(分)0510152025303540455055浓度y*10（-4）01.272.162.863.443.874.154.374.514.584.624.64解:先根据数据画图判断曲线的大致形状:观察到图形大概为一条抛物线,因此它的最高次数为2.构造其正规方程的程序如下:x0=0:5:55; x=x0; y0=0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.62 4.64; y=y0; m=2; n=length(x); b=zeros(1,m+1); f=zeros(n,m+1); for k=1:m+1f(:,k)=x.(k-1);end a=f*f; b=f*y; c=ab; c=flipud(c); disp(c); -0.0022 0.19920.2453运行上述程序得到数据的系数如下(从高次到低次排列): -0.0022 0.19920.2453对多项式系数进行重新计算,得到拟合的效果图:c=c;x=0:5:55; F=c(1)*x.2+c(2)*x+c(3);plot(x0,y0, ok,LineWidth,2),grid on;hold on;plot(x,F,-r,LineWidth,1);title(拟合后的效果图); axis(0 55 0 5)结论：掌握数据拟合的基本原理，并且掌握最小二乘法的计算方法，学会使用数学的方法对数据拟合的情况进行判断。 实验报告四题目：用种方法求数值积分摘要：按照算法计算定积分的值，要求满足给定误差的值，结果输出定积分的近视值最后的步长H前言：（目的和意义）掌握数值积分的基本思想和积分公式的推导，并写出自己设计出的变步长算法。体会事前误差和事后误差的不同之处。实验原理：Simpson公式：设计Simpson公式，即为设计含有3个节点（即复化为3阶精度）的Newton-Cotes公式。将区间划分为3等份，选取等分点作为求积节点构造求积公式，具有三阶精度的Simpson公式。将区间a,b划分为n等分，步长为h=(b-a)/n，等分点xi=a+i*h,i=0,1,2,3。设计复化求积法，先用低阶求积公式求得每个子段上的积分值，然后再将它们累加求和，用各段积分之和作为所求积分的近似值。即设计出复化Simpson公式。复化梯形公式：将区间a,b划分为n等分，步长为h=(b-a)/n，等分点xi=a+i*h,i=0,1,2,3。设计复化求积法，先用低阶求积公式求得每个子段上的积分值，再将它们累加求和，用各段积分之和作为所求积分的近似值。再根据梯形公式即可设计出复化梯形公式。Romberg公式：再简化Cotes值。将积分区间a,b分为8等份，等分点xi=a+i*(b-a)/8，i=0,1,2,8,则二分前后的Cotes值可求，再对求得的Cotes公式进行松弛，提高精度，设计出Romberg公式。三点Gauss公式：根据Newton-Cotes公式，但对求积节点自由选择，适当选取待定参数使公式具有高精度，即设计出Gauss公式。而设计出具有5阶精度即为3点Gauss公式。4用Simpson公式求积分exp(-x) dx并估计误差. a=0; b=1; f=inline(exp(-x); I=(b-a)/6*(feval(f,a)+4*feval(f,(a+b)/2)+feval(f,b)I = 0.6323336800036638.用Romberg方法计算积分2/ exp(-x) dx,要求误差不超过10(-8).编写积分法的函数M文件,源程序如下(romberg.m)function I,T=romberg(f,a,b,n,Eps)if narginEps) & (jn)| (j double(int(1/x,1,3)ans = 1.098612288668110小结：通过各种方法进行数值积分，认识各个方法的区别，通过各种题型来训练来熟练各种各种数值积分的方法。深入掌握复化Simpson公式、复化梯形公式、Romberg公式以及三点Gauss公式，知道了每种算法的关键设计要点，了解了每种算法设计的注意事项。 实验报告五题目：常微分方程数值解法摘要：科学技术中常常需要求解常微分方程的定解问题，虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法求解。前言：（目的和意义）掌握求解常微分方程的一些解法。数学原理：所谓数值解法，就是寻求解y（x）在一系列离散节点 X1x2xnx（n+1）上的近似值y1，y2，yn，y（n+1），。相邻两个节点的间距h=x（n+1）-xn称为步长。总是假定h为定数，这时节点为xn=x0+nh，n=0,1,2，。第三题：用改进的Euler方法解 取步长h=0.1计算，并与准确解y=-x-1+2相比较。解：Euler法F=x*x+x-y;a=0;b=1;h=0.1;n=(b-a)/h;X=a:h:b;Y=zeros(1,n+1);Y(1)=0;for i=2:n+1 x=X(i-1); y=Y(i-1); Y(i)=Y(i-1)+eval(F)*h;endY(6)ans =0.1186Euler校正法Y1=zeros(1,n+1);Y1(1)=0;for i=2:n+1 x=X(i-1); y=Y1(i-1); ty=Y1(i-1)+eval(F)*h;Y1(i)=Y1(i-1)+h/2*eval(F);x=X(i); y=ty; Y1(i)=Y1(i)+h/2*eval(F);endY1(6)ans =0.1450准确值temp=;f=dsolve(Dy= x*x+x-y ,y(0)=1,x);df=zeros(1,n+1);for i=1:n+1 temp=subs(f,x,X(i); df(i)=double(vpa(temp);enddisp( 步长 Euler法 Euler预测-校正公式 准确值);disp(X,Y,Y1,df)运行以上程序得到如下结果: 步长 Euler法 Euler预测-校正公式 准确值 0 0 0 1.0000 0.1000 0 0.0055 0.9100 0.2000 0.0110 0.0219 0.8400 0.3000 0.0339 0.0501 0.7900 0.4000 0.0695 0.0909 0.7600 0.5000 0.1186 0.1450 0.7500 0.6000 0.1817 0.2130 0.7600 0.7000 0.2595 0.2954 0.7900 0.8000 0.3526 0.3929 0.8400 0.9000 0.4613 0.5059 0.91001.0000 0.5862 0.6348 1.0000同时作出的对比图如下:由图看出Euler和Euler预测-校正法,都不好第六题：取h=0.2，用经典的四阶Runge-Kutta方法解下列初值问题。（1） （2）解：（1）四阶经典R-K公式作数值计算clc;F=3*y/(1+x);a=0;b=1;h=0.2;n=(b-a)/h;X=a:h:b;Y=zeros(1,n+1);Y(1)=1;for i=1:n x=X(i); y=Y(i); K1=h*eval(F); x=x+h/2; y=y+K1/2; K2=h*eval(F); x=x; y=Y(i)+K2/2; K3=h*eval(F); x=X(i)+h; y=Y(i)+K3; K4=h*eval(F); Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6;end准确解temp=;f=dsolve(Dy=3*y/(1+x),y(0)=1,x);df=zeros(1,n+1);for i=1:n+1 temp=subs(f,x,X(i); df(i)=double(vpa(temp);enddisp( 步长 四阶经典R-K法 准确值);disp(X,Y,df);figure;plot(X,df,k*,X,Y,-r);grid on;title(四阶经典R-K法解常微分方程);legend(准确值,四阶经典R-K法);运行上述程序得到如下结果 步长 四阶经典R-K法 准确值 0 1.0000 1.0000 0.2000 1.2428 1.2428 0.4000 1.5836 1.5836 0.6000 2.0442 2.0442 0.8000 2.6510 2.65111.0000 3.4365 3.4366作出的函数图形如下:（2） 四阶经典R-K公式作数值计算clc;F=3y/(1+x);a=0;b=1;h=0.2;n=(b-a)/h;X=a:h:b;Y=zeros(1,n+1);Y(1)=1;for i=1:n x=X(i); y=Y(i); K1=h*eval(F); x=x+h/2; y=y+K1/2; K2=h*eval(F); x=x; y=Y(i)+K2/2; K3=h*eval(F); x=X(i)+h; y=Y(i)+K3; K4=h*eval(F); Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6;endtemp=;f=dsolve(Dy=3y/(1+x),y(0)=1,x);df=zeros(1,n+1);for i=1:n+1 temp=subs(f,x,X(i); df(i)=double(vpa(temp);enddisp( 步长 四阶经典R-K法 准确值);disp(X,Y,df);figure;plot(X,df,k*,X,Y,-r);grid on;title(四阶经典R-K法解常微分方程);legend(准确值,四阶经典R-K法); 运行上述程序得到如下结果: 步长 四阶经典R-K法 准确值 0 1.0000 1.0000 0.2000 1.7275 1.7280 0.4000 2.7430 2.7440 0.6000 4.0942 4.0960 0.8000 5.8292 5.83201.0000 7.9960 8.0000小结：本章研究求解常微分方程的差分方法，构造数值积分主要有两条途径：基于数值积分的构造方法和基于Taylor展开的构造方法，后一种方法灵活，具有一般性，同时可得到截断误差的估计。直接用Taylor展开到处的Taylor级数法不便用于实际，给予Taylor展开构造的四阶Runge-Kutta方法则是电子计算机上的常用算法，具有精度高，程序简单，计算过程稳定，易于调节步长等优点，但也有不足之处，它所要求的函数具有教高的光滑性，且计算量大。实验报告六题目：方程求根摘要：数学中的很多问题常归结为求解数学函数方程f（x）=0，这里f（x）可以使代数多项式，也可以是超越函数。方程f（x）=0的解x*称为它的根，或者称为f（x）的零点。前言：（目的和意义）掌握方程的求根方法及原理。数学原理：设函数f（x）在a,b上连续，且f（a）f（b） format long;f=inline(1-x*sin(x);a=0;b=1;Eps=0.005;for k=1:50A(k)=a;B(k)=b;ya=feval(f,a);yb=feval(f,b);temp=(a+b)/2;X(k)=temp;yt=feval(f,temp);F(k)=yt;if abs(yt)Eps break;endif yt*ya0 a=a;b=temp;elseif yt*yb0a=temp;b=b;endend;disp(k a(k) b(k) x(k) f(x) );H=1:k,A,B,X,F;disp(H);disp(x=);disp(X(k);disp(y=);disp(yt);format shortk a(k) b(k) x(k) f(x) 1.000000000000000 0 1.000000000000000 0.500000000000000 0.760287230697899 2.000000000000000 0 1.000000000000000 0.500000000000000 0.760287230697899 3.000000000000000 0 1.000000000000000 0.500000000000000 0.760287230697899 4.000000000000000 0 1.000000000000000 0.500000000000000 0.760287230697899 5.000000000000000 0 1.000000000000000 0.500000000000000 0.760287230697899 6.000000000000000 0 1.0000000000