数值分析试题

1、1.设有某实验数据如下：0.10.20.30.40.50.60.70.80.95.12345.30535.56845.93786.42707.07987.94939.025310.3627（1）在MATLAB中作图观察离散点的结构，用最小二乘法拟合一个合适的多项式函数；（2）在MATLAB中作出拟合曲线图. 解：（1）在MATLAB中作图观察离散点的结构，用最小二乘法拟合一个合适的多项式函数。具体程序如下：>> x=0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9;>> y=5.1234 5.3053 5.5684 5.9378 6.4270 7.

2、0798 7.9493 9.0253 10.3627;>> plot(x,y,'r+'),legend('x,y'),xlabel('x'),ylabel('y'),title('(x,y)的散点结构图')从图中各点可以看到，（x,y）散点结构图的变化趋势与二次多项式很接近，故选取曲线函数为二次多项式函数，拟合多项式的次数为2。用最小二乘法拟合一个合适的多项式函数，具体程序如下：>> s=polyfit(x,y,2);>> P=poly2str(s,'t')P =

3、 8.2191 t2 - 1.8822 t + 5.3139则曲线函数为发f（x）=8.2192t2-1.8822t+5.3139（2） 在MATLAB中作出拟合曲线图，具体程序如下：>> x=linspace(0.1,0.9,9);>> y=8.2191.*x.2-1.8822.*x+5.3139;>> plot(x,y),legend('x,y'),xlabel('x'),ylabel('y'),title('拟合曲线y=f(x)')2. 在MATLAB中用复合Simpson公式编程计算下列

4、积分，并使其结果至少有6位有效数字. 解：（1）所安装的MATLAB中不含有复合辛普森程序，用M文件编辑器编写被积函数的M文件，并以jfSimpson.m命名保存至MATLAB搜索路径下。程序如下：function I,step = jfSimpson(f,a,b,type,eps)%type = 1 辛普森公式%type = 2 复合辛普森公式 if(type=2 && nargin=8) eps=1.0e-8; %缺省精度为0.0001end I=0;switch type case 1, I=(b-a)/6)*(subs(sym(f),findsym(sym(f),a)+

5、. 4*subs(sym(f),findsym(sym(f),(a+b)/2)+. subs(sym(f),findsym(sym(f),b); step=1; case 2, n=2; h=(b-a)/2; I1=0; I2=(subs(sym(f),findsym(sym(f),a)+subs(sym(f),findsym(sym(f),b)/h; while abs(I2-I1)>eps n=n+1; h=(b-a)/n; I1=I2; I2=0; for i=0:n-1 x=a+h*i; x1=x+h; I2=I2+(h/6)*(subs(sym(f),findsym(sym(f

6、),x)+. 4*subs(sym(f),findsym(sym(f),(x+x1)/2)+. subs(sym(f),findsym(sym(f),x1); end end I=I2; step=n;end （2）在命令窗口中，调用jfSimpson函数进行求解。命令窗口显示如下：>> I,s=jfSimpson(inline('exp(x.2)'),0,2,2,1.0e-8)I = 16.4526s = 102说明：求解积分函数为16.4526。3. 在MATLAB中试用经典四阶Runge-Kutta法编程求解初值问题 解：（1）用M文件编辑器将RK4.m、RK

7、_fun.m、RK_main.m的主程序命名并保存至MATLAB搜索路径下。主程序如下：RK4.mfunction t,y = RK4(func,t0,tt,y0,N,varagin)% Rk方法计算一阶级微分方程组，% 由微分方程知识可以知道，对于高阶微分方程可以化为一阶微分方程组。% 初始时刻为t0，结束时为tt，初始时刻函数值为y0% N为 步数 if nargin<4 N=100 end % 如果没有输入N的话，那么默认计算步长为（tt-t0）/100 y(1,:) = y0(:)'h = (tt-t0)/N; %步长t = t0+0:N'*h; for k =

8、1:N f1 = h*feval(func,t(k),y(k,:); f1 = f1(:)'% RK中的k1f2 = h*feval(func,t(k) + h/2,y(k,:) + f1/2); f2 = f2(:)'% RK中的k2 f3 = h*feval(func,t(k) + h/2,y(k,:) + f2/2); f3 = f3(:)'% RK中的k3f4 = h*feval(func,t(k) + h,y(k,:) + f3); f4 = f4(:)'% RK中的k4 y(k + 1,:) = y(k,:) + (f1 + 2*(f2 + f3)

9、+ f4)/6; % 所求结果endRK_fun.m%RK4方法的测试函数function f=RK_fun(t,y)f=-3*y2+2*t2+3*t;RK_main.m%rk方法的主函数 x0=0;xt=1;y0=1;%符号计算给出分析解y=dsolve('Dy=-3*y2+2*t2+3*t','y(0)=1','t') %RK4方法给出计算结果x,yrk = RK4('RK_fun',x0,xt,y0,10); %对应于真解的离散数据 yreal=subs(y,x); tol=yreal-yrk; %RK4方法与分析解的误差p

10、lot(x,yreal,x,yrk,'+',x,tol,'*')legend('分析解','RK4计算结果','两者误差') yrk;x,yrk,tol%yrk为所要计算的结果（2） 运行结果如下： >> RK_main y =<1x1 sym>ans= <1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym>

11、<1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym>

12、<1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym>4. 在MATLAB中利用Newton法编程计算方程 的一个根. 解：(1)所安装的MATLAB中不含有newton.m程序，用M文件编辑器编写被积函数的M文件，并以newtonqxf.m命名保存至MATLAB搜索路径下。程序如下：function root=newtonqxf(f,a,b,eps)% 牛顿法求函f数在区间a,b上的一个零点% f 为函数名% a 为区间左端点% b 为区间右端点% eps 为根的精度% root 为求出的函数零点

13、if(nargin=3) eps=1.0e-4;end f1=subs(sym(f),findsym(sym(f),a);f2=subs(sym(f),findsym(sym(f),b);if(f1=0) root=a;endif(f2=0) root=b;end if(f1*f2>0) disp('两端点函数值乘积大于0!'); return;else tol=1; fun=diff(sym(f); fa=subs(sym(f),findsym(sym(f),a); fb=subs(sym(f),findsym(sym(f),b); dfa=subs(sym(fun),

14、findsym(sym(fun),a); dfb=subs(sym(fun),findsym(sym(fun),b); if(dfa>dfb) root=a-fa/dfa; else root=b-fb/dfb; end while(tol>eps) r1=root; fx=subs(sym(f),findsym(sym(f),r1); dfx=subs(sym(fun),findsym(sym(fun),r1); root=r1-fx/dfx; tol=abs(root-r1); endend (2) 在命令窗口中，调用newtonqxf函数进行求解。命令窗口显示如下:>&

15、gt; root=newtonqxf('x-log(x)-2',3,1.0e-9)root =3.14625. 已知线性方程组 （1）利用MATLAB说明解该方程组的Gauss-Seidel迭代法是否收敛； （2）在MATLAB中利用Gauss-Seidel迭代法求解该方程组，并给出迭代次数. 解：通过M文件编辑器加入gauseidel.m函数程序，如下：function x,n=gauseidel(A,b,x0,eps,M)% 求解线性方程组的迭代法,其中,% A为方程组的系数矩阵;% b为方程组的右端项;% x0为迭代初始化向量% eps为精度要求,缺省值为1e-5;% M

16、为最大迭代次数,缺省值100;% x为方程组的解;% n为迭代次数;if nargin=3 eps= 1.0e-6; M = 200;elseif nargin = 4 M = 200;elseif nargin<3 error return;end D=diag(diag(A); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵G=(D-L)U;f=(D-L)b;x=G*x0+f;n=1; %迭代次数while norm(x-x0)>=eps x0=x; x=G*x0+f; n=n+1; if(n>=M) disp('Warning: 迭代次数太多，可能不收敛！'); return; endend（1）在命令窗口中，输入程序如下：>> A=5,2,1;-1,4,2;2,-5,10A = 5 2 1 -1 4 2 2 -5 10>> D=diag(diag(A)D = 5 0 0 0 4 0 0 0 10>> L=-tril(A,1)L = -5 -2 0 1 -4 -2 -2 5 -10>> U=-triu(A,1)U = 0 -2 -1 0 0 -2 0 0 0>> G=(D-L)UG = 0 -0.1945 -