《矩阵与数值分析》课程数值实验大作业.doc

《《矩阵与数值分析》课程数值实验大作业.doc》由会员分享，可在线阅读，更多相关《《矩阵与数值分析》课程数值实验大作业.doc（17页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、优质文本2011级工科硕士研究生?矩阵与数值分析?课程数值实验班 级： 学 号： 姓 名： 任课教师： 大连理工大学2011年12月20日优质文本一、 对于数列，有如下两种生成方式1、首项为，递推公式为；2、前两项为，递推公式为；给出利用上述两种递推公式生成的序列的第50项。【按第一种递推公式】clearclca=1;for i=1:50-1 a=a a(i)/3;enddisp(数列第50项小数表达为：)format longdisp(a(50)disp(分数表达为：)format ratdisp(a(50)format short运行结果数列第50项小数表达为： 4.17886670729

2、5615e-024分数表达为： 1/239299329230617530000000【按第二种递推公式】clearclca=1 1/3;for i=2:50-1 a=a 10/3*a(i)-a(i-1);endformat ratdisp(数列第50项为：)disp(a(50)format short运行结果数列第50项为：2060436 【分析】第一种算法数值稳定，计算过程舍入误差被严格控制，且按1/3的公差不断缩小。但第二种算法数值不稳定。另外，在第二种算法中，假设将递推公式“a=a 10/3*a(i)-a(i-1)中的分母移动位置，改写成“a=a 10*a(i)/3-a(i-1)，那么程

3、序运行结果为-4966040，可以舍入误差被放大的十分严重。二、 利用迭代格式及Aitken加速后的新迭代格式求方程在内的根【未经加速的代码】clceps=1e-15;i=1;x0=1;format longwhile i100 x1=sqrt(10/(x0+4); if abs(x1-x0)=eps break end x0=x1; i=i+1;enddisp(方程的解精度10(-15)disp(x1)disp(未经加速的迭代次数)disp(i)运行结果方程的解精度10(-15) 1.36523001341410未经加速的迭代次数 18【经Aitken加速的代码】clceps=1e-15;i

4、=1;x0=1;format longwhile i100 x1=sqrt(10/(x0+4); y=sqrt(10/(x1+4); z=sqrt(10/(y+4); x1=z-(z-y)2/(z-2*y+x1); if abs(x1-x0)=eps break end x0=x1; i=i+1;enddisp(方程的解精度10(-15)disp(x1)disp(未经加速的迭代次数)disp(i)运行结果方程的解精度10(-15) 1.36523001341410未经加速的迭代次数 3【分析】Aitken加速能对数列xk起明显的加速作用，在要求相同方程解精度的条件下，它能将迭代次数显著降低。实

5、际上，Aitken有时甚至能将发散的数列加速后变为收敛。三、解线性方程组 1分别Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组迭代法计算停止的条件为： 2. 用Gauss列主元消去法、QR方法求解如下方程组：【1. Jacobi方法】clci=1;eps=1e-6;A= 6 2 1 -2; 2 5 0 -2; -2 0 8 5; 1 3 2 7;b=4 7 -1 0;x0=zeros(4,1);D=diag(diag(A);L=-tril(A,-1);U=-triu(A,1);B=inv(D)*(L+U);f=inv(D)*b;while i100 x1=B*x0+f; if

6、 norm(x1-x0)=eps break end x0=x1; i=i+1;enddisp(方程的解精度10(-6)disp(x1)disp(迭代次数)disp(i)运行结果方程的解精度10(-6) 0.05204951386229 1.15094332647528 0.24463239101199 -0.57059207123432迭代次数 16【1. Gauss-Seidel方法】clci=1;eps=1e-6;A= 6 2 1 -2; 2 5 0 -2; -2 0 8 5; 1 3 2 7;b=4 7 -1 0;x0=zeros(4,1);D=diag(diag(A);L=-tril

7、(A,-1);U=-triu(A,1);B=inv(D-L)*(U);f=inv(D-L)*b;while i100 x1=B*x0+f; if norm(x1-x0)=eps break end x0=x1; i=i+1;enddisp(方程的解精度10(-6)disp(x1)disp(迭代次数)disp(i)运行结果方程的解精度10(-6) 0.05204946628111 1.15094338455986 0.24463241176981 -0.57059206335719迭代次数 8【2. Gauss列主元消去法】clcA= 2 2 1 2; 4 1 3 -1; -4 -2 0 1;

8、2 3 2 3;b=1 2 1 0;Ab=A b;N,N=size(A);x=zeros(N,1);for p=1:N-1 max1,j=max(abs(A(p:N,p); temp=Ab(p,:); Ab(p,:)=Ab(j+p-1,:); Ab(j+p-1,:)=temp; if Ab(p,p)=0 disp(方程无解) break end for k=p+1:N mult=Ab(k,p)/Ab(p,p); Ab(k,p:N+1)=Ab(k,p:N+1)-mult*Ab(p,p:N+1); endendb=Ab(:,N+1);xx=zeros(N,1);for k=N:-1:1 xx(k)

9、=b(k)/Ab(k,k); i=(1:k-1); b(i)=b(i)-xx(k)*Ab(i,k);endx=xx运行结果x = 1.54166666666667 -2.75000000000000 0.08333333333333 1.66666666666667【2. QR分解法】clcfor i=1:3 Ai=zeros(5-i); Qi=eye(4);endx=zeros(4,1);y=zeros(4,1);R=zeros(4);A1=2 2 1 2; 4 1 3 -1; -4 -2 0 1; 2 3 2 3;b=1,2,1,0;for i=1:3 I=eye(size(Ai); a=

10、Ai(:,1); w=a-norm(a)*I(:,1); hw=I-2/(w*w)*(w*w); Qi(i:4,i:4)=hw; if i=3 break end c=hw*Ai; Ai+1=c(2:max(size(Ai),2:max(size(Ai);endQz=(Q3*Q2*Q1);R=Q3*Q2*Q1*A1;c=Qz*b;x(4)=c(4)/R(4,4);for i=3:-1:1 x(i)=(c(i)-R(i,i+1:4)*x(i+1:4)/R(i,i);endx运行结果x = 1.54166666666666 -2.75000000000000 0.08333333333333 1.

11、66666666666666【分析】Jacobi迭代法和Gauss-Seidel迭代法都可用来求解线性方程组。在同等精度下，求解本道题Jacobi迭代法迭代了16次，而Gauss-Seidel仅为8次，是前者的一半。所以可以认为，在某些情况下，Gauss-Seidel是比拟好的解法，但不总如此。Gauss消去法可能发生小主元做除数，从而导致计算结果严重偏离真实值，所以在消元过程中，每一步都按列选主元，也就是Gauss列主元消去法，它可以有效防止过于严重的舍入误差。正交矩阵是性态最好的矩阵，它不改变矩阵的条件数，通过QR分解来计算线性方程组，也可以防止误差的放大，保证计算过程具有数值稳定性。四、

12、一组数据点，编写一程序求解三次样条插值函数满足 并针对下面一组具体实验数据0.250.30.390.450.530.50000.54770.62450.67080.7280求解，其中边界条件为.【程序代码，文件名 Spline.m】function s=Spline(x,y,f0,fn) % 输入实验数据x,y；边界二阶导数f0,fnclcfigure(1)plot(x,y,*r)hold onN=max(size(x);syms t s;for k=1:N-1; h(k)=x(k+1)-x(k);endg(1)=3*(y(2)-y(1)/h(1)-h(1)/2*f0;g(N)=3*(y(N)

13、-y(N-1)/h(N-1)+h(N-1)/2*fn;for k=2:N-1 lamda(k)=h(k)/(h(k)+h(k-1); miu(k)=h(k-1)/(h(k)+h(k-1);g(k)=3*(miu(k)*(y(k+1)-y(k)/h(k)+.lamda(k)*(y(k)-y(k-1)/h(k-1);endc=1,miu(2:N-1);b=2*ones(1,N);a=lamda(2:N-1),1;A=diag(c,1)+diag(b,0)+diag(a,-1);d=c;a=0,a;u(1)=b(1);for i=2:N l(i)=a(i)/u(i-1); u(i)=b(i)-l(i

14、)*c(i-1);endL=eye(N)+diag(l(2:N),-1);U=diag(u)+diag(d,1);z(1)=g(1);for i=2:N z(i)=g(i)-l(i)*z(i-1);endm(N)=z(N)/u(N);for i=N-1:-1:1 m(i)=(z(i)-c(i)*m(i+1)/u(i);endfor i=1:N-1 s(i)=(h(i)+2*(t-x(i)*(t-x(i+1)2*y(i)/(h(i)3+.(h(i)-2*(t-x(i+1)*(t-x(i)2*y(i+1)/(h(i)3+.(t-x(i)*(t-x(i+1)2*m(i)/(h(i)2+.(t-x(i

15、+1)*(t-x(i)2*m(i+1)/(h(i)2;Enddigits(8)s=vpa(expand(s) % 输出分段表达式for i=1:N-1 % 绘图 v=x(i):0.005:x(i+1); y=subs(s(i),v); plot(v,y); hold on;endgrid on命令窗口输入x=0.25 0.3 0.39 0.45 0.53;y=0.5000,0.5477,0.6245,0.6708,0.7280;f0=0;fn=0;Spline(x,y,f0,fn)运行结果ans =4.6988737*t2-.20505552*t+.35547747-6.2651650*t3,

16、-2.6329843*t2+1.9945019*t+.13552173+1.8813439*t3,.10638708*t2+.92614706*t+.27440786-.45999912*t3,-3.4093028*t2+2.5082075*t+.37098798e-1+2.1442156*t3【分析】运行结果是一个四个元素的矩阵，各个元素依次代表四个分段区间上的三次样条曲线。例如在第一个区间0.25 0.3上，拟合得到的三次样条曲线方程4.6988737*t2-0.20505552*t+0.35547747-6.2651650*t3。由图像可知，三次样条曲线是很光滑的曲线，拟合效果很好。五、

17、编写程序构造区间上的以等分结点为插值结点的Newton插值公式，假设结点数为包括两个端点，给定相应的函数值，插值区间和等分的份数，该程序能快速计算出相应的插值公式。以，为例计算其对应的插值公式，分别取不同的值并画出原函数的图像以及插值函数的图像，观察当增大时的逼近效果.【程序代码，文件名 Newton.m】function f1=Newton(n)clcsyms x t fa=-1;b=1;f=1./(1+25*x.2);y=zeros(1,n);x=linspace(a,b,n); h=-1:0.01:1;m=n;c(1:m) = 0.0;yh=subs(f,h);y=subs(f,x);

18、yk=y;f1=y(1); y1=0;k=1;for i=1:m-1 for j=i+1:m y1(j)=(y(j)-y(i)/(x(j)-x(i); end c(i)=y1(i+1); k=k*(t-x(i); f1=f1+c(i)*k; y=y1;endf1=simplify(f1) yfh=subs(f1,h);figure(1)plot(h,yfh,-k,x,yk,*r)grid onhold onx=-1:0.01:1;y=1./(1+25*x.2);plot(x,yh,b)legend(插值曲线,插值节点,被插曲线);运行结果Newton(4) 5个插值节点拟合方程为：1-3225/754*t2+1250/377*t4Newton(5) 6个插值节点拟合方程为： 用digits(4),vpa(ans)处理后.5673+.1599e-16*t-1.731*t2-.4441e-15*t3+1.202*t4+.1110e-14*t5当n继续增大时有：n=7n=10n=12n=15【分析】此题说明在不分段拟合的条件下，节点并不是越不越好。插值节点增多会导致插值多项式的次数增高，而高次多项式的振荡次数增多有可能使插值多项式在非节点处的误差变得很大Runge现象。所以在高次插值时，为了克服这种现象，建议改为分段低次插值。