数值分析上机实验报告.doc

西安工程大学数值计算方法实验报告 实验报告一题目： （绪论） 非线性方程求解及误差估计摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法、Newton法和改进的Newton法。可以节省计算机的计算时间，还能减小不必要的误差。前言：（目的和意义）掌握二分法与Newton法的基本原理、应用以及熟练掌握用MATLAB求函数积分数学原理：（1） 函数的调用格式：quadl(filename,a,b,tol,trace) 其中filename是调用函数名，a和b分别为定积分的下限和上限。用来控制积分精度。（2）秦九韶算法： Sn=an Sk=xSk+1+ak (k=n-1,n-2,.,0), Pn(x)=S0程序设计：例1.1 计算积分 利用MATLAB，下面给出主程序g=inline(x.10.*exp(x-1); %定义一个语句函数g(x)=exp(x10*exp(x-1)I=quadl(g,0,1) I =0.0098例1.9 秦九韶算法a0=3,ak=2ak-1+3,Pn(x)=anxn+an-1x(n-1)+.+a1x+a0求I1=P100(0.5),I2=P150(13)x=input(x=); n=input(n=); a=3;for i=1:n a=2*a+3;end s=z; b=(a-3)/2; for m=1:100s=x*s+b; b=(b-3)/2; end disp(s);x=0.5n=100 600.0000x=3n=100 4.7039e+078结果分析和讨论：结论：对于二分法，只要能够保证在给定的区间内有根，使能够收敛的，当时收敛的速度和给定的区间有关，二且总体上来说速度比较慢。Newton法，收敛速度要比二分法快，但是最终其收敛的结果与初值的选取有关，初值不同，收敛的结果也可能不一样，也就是结果可能不时预期需要得结果。改进的Newton法求解重根问题时，如果初值不当，可能会不收敛，这一点非常重要，当然初值合适，相同情况下其速度要比Newton法快得多。实验报告二题目： (插值法)用各类插值方法法求解相关数值分析问题摘要： 非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。前言：（目的和意义）掌握插值法的基本原理及用插值法求解相关数值分析问题的方法以及熟练掌握用MATLAB求函数积分。数学原理： MATLAB提供了一、三、N维数据插值函数以及三次样条插值函数等。数据插值问题为一维插值，其采用的方法有线性方法、最近方法、三次样条和三次插值，在MATLAB中，实现这些插值的函数为interp1，其调用格式为inpert1(a,b,c,method)；函数根据a、b的值，计算函数在c处的值，method为插值方法，linear 为线性插值，cubic为三次多项式插值，spline 为三次样条插值。程序设计：3.给出f(x)=lnx,用一、二、三次线形插值计算ln0.54的近似值x0.40.50.60.70.8-0.916291-0.693147-0.510826-0.357765-0.223144x=0.4:0.1:0.8;f=-0.916291,-0.693147,-0.510826,-0.357765,-0.223144;I1=interpl(x,f,0.54) %用一次线形插值计算f(x) I1 = -0.620218600000000I3=interp1(x,f,0.54,spline) %用次样条插值计算f(x)I3 = -0.61597777000000021设f(x)=1/(1+x2),在-5，5上取n=10，安等距节点求分段线性插值函数Ib（x），计算节点间中点处的 与f（x）的值，并估计误差。 x=linspace(-5,5,10); y=1(1+(x.2); p=polyfit(x,y,1)p = 0.0000 11.18524.给定数据如下表，试求三次样条插值及Hermite插值S(x)，并满足条件: (1)S(0.25)=1.0000, S(0.53)=0.6868 (2)S(0.25)=S(0.53)=0xi0.250.300.390.450.53yi0.500000.54770.62450.67080.7280（1）x=0.25,0.30,0.39,0.45,0.53; y=0.5000,0.5477,0.6245,0.6708,0.7280; p=polyfit(x,y,spline) %三次样条插值函数程序p =0.732659297655720 -1.394857855010107 1.553581221803837 0.187348425537683（2）用Hermite程序： n=length(x0); m=length(x); for k=1:m y=0.0; for i=1:n h=1.0; a=0.0; for j=1:n if j=i h=h*(x(k)-x0(j)/(x0(i)-x0(j)2; a=1/(x0(i)-x0(j)+a; end end y=y+h*(x0(i)-x(k)*(2*a*y0(i)-dy(i)+y0(i); end y(k)=y;endx0=0.25,0.30,0.39,0.45,0.53;Y0=0.5000,0.5477,0.6245,0.6708,0.7280;H=hermite(x0,y0)结果分析和讨论：结论：实验报告三题目： （函数逼近与计算）用曲线拟合的最小二乘法求解相关数值分析问题摘要：与数值差值类似，曲线拟合的目的也是用一个简单的函数去逼近一个复杂的或未知的函数，在构造函数去逼近，不要求曲线严格通过采样点，但希望能尽量的靠近这些点，从而使误差尽可能的小。前言：（目的和意义）掌握曲线拟合的最小二乘法的基本原理以及熟练掌握用MATLAB求曲线拟合求解相关数值分析问题的方法。数学原理：程序设计：(1)观察物体的直线运动，得出时间t与距离s的关系如下表3.7，求运动方程。表3.7t/s00.91.93.03.95.0S/m010305080110t=0,0.9,1.9,3.0,3.9,5.0;s=0,10,30,50,80,110;plot(t,s) %先用MATLAB绘出图形 可以看出s是关于t的近似一元二次函数，用MATLAB曲线拟合运动方程：t=0,0.9,1.9,3.0,3.9,5.0;s=0,10,30,50,80,110;P=polyfit(x,y,2)P=2.2488,11.0814,-0.5834即运动方程为：s=2.2488x2+11.0814x-0.5834 (2)在某化学反应里，根据实验所得分解物的质量分数y于是间t的关系如表3.8所示，用最小二乘拟合求y=F(x).t/min0510152025y/01.272.162.863.443.87t/min303540455055y/4.154.374.514.584.624.64t=0:5:55; y=0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.62,4.64y = 0 1.2700 2.1600 2.8600 3.4400 3.8700 4.1500 4.3700 4.5100 4.5800 4.6200 4.6400 plot(t,y) %运用MATLAB绘出y关于t的图形为一元二次函数 P=polyfit(t,y,3) %三次、二次、一次和常数项的系数P = 0.0000 -0.0049 0.2557 0.0442y=F(t)=-0.0049+0.2557t+0.0442结论： 曲线拟合最小二乘法是数值分析研究中一个重要领域之一，通过借助于其他的数学工具，我们可以更加深刻的理解“拟合”的真实含义，求解数学问题更加容易。实验报告四题目： （数值积分与数值微分）用Simpson公式,Romberg方法,三点及五点Gauss公式,复化两点Gauss公式计算某些函数的定积分。摘要：计算某些函数的积分时，由于寻找原函数比较困难，需结合MATLAB就可以很快而且容易的算出结果。前言：（目的和意义）熟悉Simpson公式,Romberg方法,三点及五点Gauss公式,复化两点Gauss公式的形式以及熟练掌握并结合MATLAB计算某些函数。数学原理：程序设计：4.用 simpson 公式求积分 并估计误差 a=0; b=1; h=inline(exp(-x); s=(b-a)/6)*h(a)+4*h(a+b)/2)+h(b)s =0.6323 %用simpson 公式求积分 f=inline(exp(-x); S=quadl(g,0,1)S = 0.6321 %输出结果准确值误差：s-s=0.0002=10-48.用romberg 方法计算积分 ，要求误差不超过10-8if narginEps) & (jn)| (jf=inline(1/(t+2);a=-1; %求积区间b=1; f=inline(1/(t+2)f = Inline function: f(t) = 1/(t+2) I3=(5/9)*f(-sqrt(15)/5)+(8/9)*f(0)+(5/9)*(sqrt(15)/5)I3 = 1.3281 I5=0.2369268*(1/(2-0.9061798)+1/(2+0.9061798)+0.4786289*(1/(2-0.5384693)+1/(2+0.5384693)+0.5688889*1/2I5 = 1.0986（3）取m=2,n=1,2应用复化公式; 取m=2,n=3,4应用复化公式通过查找各种求积的求积节点和求积系数,按各积分公式进行计算：a=1;b=3;f=inline(1/x); %m=2,n=1Gauss-Legendre复化求积公式n=1;m=2;h=(b-a)/(m); %步长x=0.5773502692 -0.5773502692; %n=1的Gauss-Legendre求积节点;A=1 1; %n=1的Gauss-Legendre求积系数k=1;I(k)=0;for j=1:n+1 %按Gauss-Legendre公式计算 temp=0; for i=1:m tx=a+(2*(i-1)+1)*h/2+h*x(j)/2; temp=temp+feval(f,tx); end I(k)=I(k)+A(j)*temp;endI(k)=I(k)*h/2; %m=2,n=2Gauss-Legendre复化求积公式k=k+1;n=2;m=2;h=(b-a)/(m);x=0.7745966692 -0.7745966692 0; %n=2;%Gauss-Legendre求积节点A=0.5555555556 0.5555555556 0.888888889; %Gauss-Legendre求积系数I(k)=0;for j=1:n+1 %按Gauss-Legendre公式计算temp=0;for i=1:m tx=a+(2*(i-1)+1)*h/2+h*x(j)/2; temp=temp+feval(f,tx);endI(k)=I(k)+A(j)*temp;endI(k)=I(k)*h/2; %m=2,n=3 Gauss-Lobatto公式k=k+1;m=2;n=3;h=(b-a)/(m); %步长x=1 -1 0.447214 -0.447214; %求积节点A=1/6 1/6 5/6 5/6; %求积系数I(k)=0;L1=0;for i=2:m %按Gauss-Lobatto公式进行计算tx=a+(i-1)*h;L1=L1+feval(f,tx);endL1=L1*2;L1=L1+feval(f,a)+feval(f,b);L1=L1*A(1);L2=0;for j=2:n %按Gauss-Lobatto公式计算temp=0;for i=1:m tx=a+(2*(i-1)+1)*h/2+h*x(j+1)/2; temp=temp+feval(f,tx);endL2=L2+A(j+1)*temp;endI(k)=(L1+L2)*h/2; %得到积分结果k=k+1;m=2;n=4;h=(b-a)/(m); %步长x=1 -1 0 0.654654 -0.654654 ; %节点A=1/10 1/10 32/45 49/90 49/90; %系数I(k)=0;L1=0; %进入计算过程,按Gauss-Lobatto公式计算for i=2:mtx=a+(i-1)*h;L1=L1+feval(f,tx);endL1=L1*2;L1=L1+feval(f,a)+feval(f,b);L1=L1*A(1);L2=0;for j=2:ntemp=0;for i=1:m tx=a+(2*(i-1)+1)*h/2+h*x(j+1)/2; temp=temp+feval(f,tx);endL2=L2+A(j+1)*temp;endI(k)=(L1+L2)*h/2; %得到积分结果I %显示各种求积公式得到的结果 I = 1.0977 1.0986 1.0986 1.0986准确值为: double(int(1/x,1,3)ans = 1.0986结果分析和讨论：结论：数值积分可以归结为函数值的四则运算，从而使计算过程可以在计算机上完成，处理数值积分的基本方法就是逼近法，计算中使用的公式算法都是依据这个思想结合插值原理得来的。通过在计算机上编程，简化了计算，也是结果更为准确。实验报告五题目： （常微分数值方程解法）摘要：前言：（目的和意义）熟悉Euler方法解题以及熟练掌握MATLAB。数学原理：程序设计：3.用改进的Euler方法解y=x2+x-y；y(0)=0取步长h=0.1计算y(0.5)，并计算准确值y=-exp(-x)+x2-x+1. F=x2+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+1x=X(i-1);y=Y(i-1);Y(i)=Y(i-1)+eval(F)*h;endYY = 0 0 0.0110 0.0339 0.0695 0.1186 0.1817 0.2595 0.3526 0.4613 0.5862Y1=zeros(1,n+1); %Euler预测、校正公式Y1(1)=0;for i=2:n+1x=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);endYY = 0 0 0.0110 0.0339 0.0695 0.1186 0.1817 0.2595 0.3526 0.4613 0.5862temp=;f=dsolve(Dy=x2+x-y ,y(0)=0,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);figure;plot(X,df,k-,X,Y,-r,X,Y1,.-b); %通过画图来观察效果grid on;title(Euler法和Euler预测、校正法解常微分方程); legend(准确值,Euler法,Euler预测-校正法); 步长 Euler法 Euler预测-校正公式 准确值 0 0 0 0- 36 -.0.1000 0 0.0055 0.0052 0.2000 0.0110 0.0219 0.0213 0.3000 0.0339 0.0501 0.0492 0.4000 0.0695 0.0909 0.0897 0.5000 0.1186 0.1450 0.1435 0.6000 0.1817 0.2130 0.2112 0.7000 0.2595 0.2954 0.2934 0.8000 0.3526 0.3929 0.3907 0.9000 0.4613 0.5059 0.5034 1.0000 0.5862 0.6348 0.63216.取h=0.2，用经典的四阶Runge-Kutta方法求解下列方程组初值问题。（1）y=x+y,0x1; y(0)=1 (2)y=3y/(1+x),0x1; y(0)=1 (1).F=x+y; %四阶经典R-K公式作数值计算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=x2+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( 步长 四阶经典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 0.8400 0.4000 1.5836 0.7600 0.6000 2.0442 0.7600 0.8000 2.6510 0.8400 1.0000 3.4365 1.0000 (2). 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;endtemp=;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.7275 1.7280 0.4000 2.7430 2.7440 0.6000 4.0942 4.0960 0.8000 5.8292 5.8320 1.0000 7.9960 8.0000结果分析和讨论：结论：数值积分可以归结为函数值的四则运算，从而使计算过程可以在计算机上完成，处理数值积分的基本方法就是逼近法，计算中使用的公式算法都是依据这个思想结合插值原理得来的。通过在计算机上编程，简化了计算，也是结果更为准确。实验报告六题目： （方程求根）前言：（目的和意义）熟练比例求根法以及熟练掌握用MATLAB求函数积分。数学原理：程序设计：2.用比例求根法求f(x)=1-sinx=0 在区间（0，1）内的一个根，直到近似根X满足精度f(X)0.005时终止计算。format short;f=inline(1-x*sin(x);a=0;b=1;Eps=1E-5;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 short k a(k) b(k) x(k) f(x) 1.0000 0 1.0000 0.5000 0.7603 2.0000 0 1.0000 0.5000 0.7603 3.0000 0 1.0000 0.5000 0.7603 4.0000 0 1.0000 0.5000 0.7603 5.0000 0 1.0000 0.5000 0.7603 6.0000 0 1.0000 0.5000 0.7603 7.0000 0 1.0000 0.5000 0.7603 8.0000 0 1.0000 0.5000 0.7603 9.0000 0 1.0000 0.5000 0.7603 10.0000 0 1.0000 0.5000 0.7603 11.0000 0 1.0000 0.5000 0.7603 12.0000 0 1.0000 0.5000 0.7603 13.0000 0 1.0000 0.5000 0.7603 14.0000 0 1.0000 0.5000 0.7603 15.0000 0 1.0000 0.5000 0.7603 16.0000 0 1.0000 0.5000 0.7603 17.0000 0 1.0000 0.5000 0.7603 18.0000 0 1.0000 0.5000 0.7603 19.0000 0 1.0000 0.5000 0.7603 20.0000 0 1.0000 0.5000 0.7603 21.0000 0 1.0000 0.5000 0.7603 22.0000 0 1.0000 0.5000 0.7603 23.0000 0 1.0000 0.5000 0.7603 24.0000 0 1.0000 0.5000 0.7603 25.0000 0 1.0000 0.5000 0.7603 26.0000 0 1.0000 0.5000 0.7603 27.0000 0 1.0000 0.5000 0.7603 28.0000 0 1.0000 0.5000 0.7603 29.0000 0 1.0000 0.5000 0.7603 30.0000 0 1.0000 0.5000 0.7603 31.0000 0 1.0000 0.5000 0.7603 32.0000 0 1.0000 0.5000 0.7603 33.0000 0 1.0000 0.5000 0.7603 34.0000 0 1.0000 0.5000 0.7603 35.0000 0 1.0000 0.5000 0.7603 36.0000 0 1.0000 0.5000 0.7603 37.0000 0 1.0000 0.5000 0.7603 38.0000 0 1.0000 0.5000 0.7603 39.0000 0 1.0000 0.5000 0.7603 40.0000 0 1.0000 0.5000 0.7603 41.0000 0 1.0000 0.5000 0.7603 42.0000 0 1.0000 0.5000 0.7603 43.0000 0 1.0000 0.5000 0.7603 44.0000 0 1.0000 0.5000 0.7603 45.0000 0 1.0000 0.5000 0.7603 46.0000 0 1.0000 0.5000 0.7603 47.0000 0 1.0000 0.5000 0.7603 48.0000 0 1.0000 0.5000 0.7603 49.0000 0 1.0000 0.5000 0.7603 50.0000 0 1.0000 0.5000 0.7603x= 0.5000y= 0.76034.比较以下两种求exp(x)+10x-2=0的根到三位小数所需的计算量。（1）.在区间（0，1）内用二分法；（2）.用迭代法，取x0=0(1).format short;f=inline(exp(x)+10*x-2);a=0;b=1;Eps=1E-5;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=te