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

2011级工科硕士研究生《矩阵与数值分析》课程数值实验班级：学号：姓名：任课教师：大连理工大学2011年12月20日1一、对于数列11111,,,,,392781，有如下两种生成方式1、首项为01a，递推公式为11,1,2,3nnaan；2、前两项为0111,3aa，递推公式为1210,2,3,3nnnaaan；给出利用上述两种递推公式生成的序列的第50项。【按第一种递推公式】clearclca=1;fori=1:50-1a=[aa(i)/3];enddisp('数列第50项小数表达为：')formatlongdisp(a(50))disp('分数表达为：')formatratdisp(a(50))formatshort[运行结果]数列第50项小数表达为：4.178866707295615e-024分数表达为：1/239299329230617530000000【按第二种递推公式】clearclca=[11/3];fori=2:50-12a=[a10/3*a(i)-a(i-1)];endformatratdisp('数列第50项为：')disp(a(50))formatshort[运行结果]数列第50项为：2060436【分析】第一种算法数值稳定，计算过程舍入误差被严格控制，且按1/3的公差不断缩小。但第二种算法数值不稳定。另外，在第二种算法中，若将递推公式“a=[a10/3*a(i)-a(i-1)]”中的分母移动位置，改写成“a=[a10*a(i)/3-a(i-1)]”，则程序运行结果为-4966040，可以舍入误差被放大的十分严重。二、利用迭代格式110,0,1,2,4kkxkx及Aitken加速后的新迭代格式求方程324100xx在[1,1.5]内的根【未经加速的代码】clceps=1e-15;i=1;x0=1;formatlongwhilei100x1=sqrt(10/(x0+4));ifabs(x1-x0)=epsbreakend3x0=x1;i=i+1;enddisp('方程的解[精度10^(-15)]')disp(x1)disp('未经加速的迭代次数')disp(i)[运行结果]方程的解[精度10^(-15)]1.36523001341410未经加速的迭代次数18【经Aitken加速的代码】clceps=1e-15;i=1;x0=1;formatlongwhilei100x1=sqrt(10/(x0+4));y=sqrt(10/(x1+4));z=sqrt(10/(y+4));x1=z-(z-y)^2/(z-2*y+x1);ifabs(x1-x0)=epsbreakendx0=x1;i=i+1;enddisp('方程的解[精度10^(-15)]')disp(x1)disp('未经加速的迭代次数')disp(i)4[运行结果]方程的解[精度10^(-15)]1.36523001341410未经加速的迭代次数3【分析】Aitken加速能对数列｛xk｝起明显的加速作用，在要求相同方程解精度的条件下，它能将迭代次数显著降低。实际上，Aitken有时甚至能将发散的数列加速后变为收敛。三、解线性方程组1．分别Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组12346212425027,2085113270xxxx迭代法计算停止的条件为：6)()1(3110maxkjkjjxx．2.用Gauss列主元消去法、QR方法求解如下方程组：12342212141312.4201123230xxxx【1.Jacobi方法】clci=1;eps=1e-6;A=[621-2;250-2;5-2085;1327];b=[47-10]';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;whilei100x1=B*x0+f;ifnorm(x1-x0)=epsbreakendx0=x1;i=i+1;enddisp('方程的解[精度10^(-6)]')disp(x1)disp('迭代次数')disp(i)[运行结果]方程的解[精度10^(-6)]0.052049513862291.150943326475280.24463239101199-0.57059207123432迭代次数16【1.Gauss-Seidel方法】6clci=1;eps=1e-6;A=[621-2;250-2;-2085;1327];b=[47-10]';x0=zeros(4,1);D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);B=inv(D-L)*(U);f=inv(D-L)*b;whilei100x1=B*x0+f;ifnorm(x1-x0)=epsbreakendx0=x1;i=i+1;enddisp('方程的解[精度10^(-6)]')disp(x1)disp('迭代次数')disp(i)[运行结果]方程的解[精度10^(-6)]0.052049466281111.150943384559860.24463241176981-0.570592063357197迭代次数8【2.Gauss列主元消去法】clcA=[2212;413-1;-4-201;2323];b=[1210]';Ab=[Ab];[N,N]=size(A);x=zeros(N,1);forp=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;ifAb(p,p)==0disp('方程无解')breakendfork=p+1:Nmult=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);fork=N:-1:18xx(k)=b(k)/Ab(k,k);i=(1:k-1)';b(i)=b(i)-xx(k)*Ab(i,k);endx=xx[运行结果]x=1.54166666666667-2.750000000000000.083333333333331.66666666666667【2.QR分解法】clcfori=1:3A{i}=zeros(5-i);Q{i}=eye(4);endx=zeros(4,1);y=zeros(4,1);R=zeros(4);A{1}=[2212;413-1;-4-201;2323;];b=[1,2,1,0]';fori=1:3I=eye(size(A{i}));a=A{i}(:,1);w=a-norm(a)*I(:,1);hw=I-2/(w'*w)*(w*w');Q{i}(i:4,i:4)=hw;ifi==39breakendc=hw*A{i};A{i+1}=c(2:max(size(A{i})),2:max(size(A{i})));endQz=(Q{3}*Q{2}*Q{1})';R=Q{3}*Q{2}*Q{1}*A{1};c=Qz'*b;x(4)=c(4)/R(4,4);fori=3:-1:1x(i)=(c(i)-R(i,i+1:4)*x(i+1:4))/R(i,i);endx[运行结果]x=1.54166666666666-2.750000000000000.083333333333331.66666666666666【分析】Jacobi迭代法和Gauss-Seidel迭代法都可用来求解线性方程组。在同等精度下，求解本道题Jacobi迭代法迭代了16次，而Gauss-Seidel仅为8次，是前者的一半。所以可以认为，在某些情况下，Gauss-Seidel是比较好的解法，但不总如此。Gauss消去法可能发生小主元做除数，从而导致计算结果严重偏离真实值，所以在消元过程中，每一步都按列选主元，也就是Gauss列主元消去法，它可以有效避免过于严重的舍入误差。正交矩阵是性态最好的矩阵，它不改变矩阵的条件数，通过QR分解来计算线性方程组，也可以避免误差的放大，保证计算过程具有数值稳定性。四、已知一组数据点，编写一程序求解三10次样条插值函数满足并针对下面一组具体实验数据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')holdonN=max(size(x));symsts;fork=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)-y(N-1))/h(N-1)+h(N-1)/2*fn;fork=2:N-1lamda(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);11d=c;a=[0,a];u(1)=b(1);fori=2:Nl(i)=a(i)/u(i-1);u(i)=b(i)-l(i)*c(i-1);endL=eye(N)+diag(l(2:N),-1);U=diag(u)+diag(d,1);z(1)=g(1);fori=2:Nz(i)=g(i)-l(i)*z(i-1);endm(N)=z(N)/u(N);fori=N-1:-1:1m(i)=(z(i)-c(i)*m(i+1))/u(i);endfori=1:N-1s(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+1))*(t-x(i))^2*m(i+1)/(h(i))^2;Enddigits(8)s=vpa(expand(s))%输出分段表达式fori=1:N-1%绘图v=x(i):0.005:x(i+1);y=subs(s(i),v);plot(v,y);holdon;endgridon12[命令窗口输入]x=[0.250.30.390.450.53];y=[0.5000,0.5477,0.6245,0.6708,0.7280];f0=0;fn=0;Spline(x,y,f0,fn)[运行结果]ans=[4.6988737*t^2-.20505552*t+.35547747-6.2651650*t^3,-2.6329843*t^2+1.9945019*t+.13552173+1.8813439*t^3,.10638708*t^2+.92614706*t+.27440786-.45999912*t^3,-3.4093028*t^2+2.5082075*t+.37098798e-1+2.1442156*t^3]【分析】运行结果是一个四个元素的矩阵，各个元素依次代表四个分段区间上的三次样条曲线。例如在第一个区间[0.250.3]上，拟合得到的三次样条