大工矩阵编程实验报告.docx

2011级工科硕士研究生矩阵与数值分析课程数值实验2011级工科硕士研究生矩阵与数值分析课程数值实验一、 对于数列，有如下两种生成方式1、首项为，递推公式为；2、前两项为，递推公式为；给出利用上述两种递推公式生成的序列的第50项。第一种方法：clc;clear;a=1;for i=2:50 a=1/3*a; i=i+1;enddisp(第50项为:)disp(a)结果如下：第50项为: 4.1789e-024第二种方法：clc;clear;a=1;b=1/3;for i=3:50 c=10/3*b-a; a=b; b=c; i=i+1;enddisp(第50项为:)disp(c)结果如下：第50项为: 2.0604e+006分析：由上述结果可知，第一种递推公式计算结果是精确的，而第二种递推公式计算的结果前几项比较接近，后来计算的误差就越来较大了。二、 利用迭代格式及Aitken加速后的新迭代格式求方程在内的根Aitken加速后：clc;clear;format long a=1;b=sqrt(10/(a+4);c=sqrt(10/(b+4);x=c-(c-b)2/(a-2*b+c);i=0;while(x3+4*x2-10=0) d=sqrt(10/(c+4); a=b; b=c; c=d; x=c-(c-b)2/(a-2*b+c); i=i+1;enddisp(此方程的根为：)disp(x)disp(迭代次数为：)disp(i)结果如下：此方程的根为： 1.36523001341410迭代次数为： 7未加速：clc;clear;format long x=1; i=0;while(x3+4*x2-10=0) x=sqrt(10/(x+4); i=i+1;enddisp(此方程的根为：)disp(x)disp(迭代次数为：)disp(i)结果如下：此方程的根为： 1.36523001341410迭代次数为： 18三、解线性方程组 1分别Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组迭代法计算停止的条件为：disp(此线性方程组的解为：x=)disp(y)disp(迭代次数为：)disp(k)结果如下：此线性方程组的解为：x= 0.05204976365435 1.15094341803680 0.24463224043774 -0.57059192015745迭代次数为： 15Jacobi迭代法：clc;clear;format long A=6 2 1 -2;2 5 0 -2; -2 0 8 5;1 3 2 7;b=4 7 -1 0;delta=10(-6);%误差n=length(b);k=1;x=zeros(n,1);y=zeros(n,1);while 1 for i=1:n y(i)=b(i); for j=1:n y(i)=y(i)-A(i,j)*x(j); end y(i)=x(i)+y(i)/A(i,i); end if norm(y-x,inf)delta break; end x=y; k=k+1;end Gauss-Seidel迭代法：clc;clear;format long A=6 2 1 -2;2 5 0 -2; -2 0 8 5;1 3 2 7;b=4 7 -1 0;delta=10(-6);%误差n=length(b);k=1;x=zeros(n,1);y=zeros(n,1);while 1 for i=1:n y(i)=b(i); for j=1:i-1 y(i)=y(i)-A(i,j)*y(j); end for j=i:n y(i)=y(i)-A(i,j)*x(j); end y(i)=x(i)+y(i)/A(i,i); end if norm(y-x,inf)B(max,k) max=i; end end temp=B(k,k:n+1); B(k,k:n+1)=B(max,k:n+1); B(max,k:n+1)=temp; for i=k+1:n B(i,k)=-B(i,k)/B(k,k); B(i,k+1:n+1)=B(i,k+1:n+1)+B(i,k)*B(k,k+1:n+1); endendx(n,1)=B(n,n+1)/B(n,n); QR方法：clc;clear;a=2 2 1 2;4 1 3 -1; -4 -2 0 1;2 3 2 3;b=1;2;1;0;x=qrfenjie(a,b);disp(此方程组的解为：x=)disp(x)单独m文件：function x=qrfenjie(b,c)h=b;q=zeros(max(size(b),max(size(b)*(max(size(b)-1);g=max(size(b);for i=1:max(size(b)-1 z=eye(size(b); w=b(:,1)-norm(b(:,1)*z(:,1); q(1:max(size(b),1+(i-1)*g:(i-1)*g+max(size(b)=z-2/(w*w)*w*w; b=q(1:max(size(b),1+(i-1)*g:(i-1)*g+max(size(b)*b; b=b(2:max(size(b),2:max(size(b);endfor i=1:g-2q(1:g,1:g)=eye(i,i),zeros(i,g-i);zeros(g-i,i),q(1:g-i,i*g+1:i*g+g-i)*q(1:g,1:g);endd=q(1:g,1:g)*h,c;x=zeros(1,g);a1=d(:,1:g);a2=d(:,g+1);x(g)=a2(g)/a1(g,g);for k = g-1:-1:1 x(k)=a2(k); for p=g:-1:k+1 x(k)= x(k)-a1(k,p)*x(p); end x(k)=x(k)/a1(k,k); end结果如下：此方程组的解为：x= 1.5417 -2.7500 0.0833 1.6667四、已知一组数据点，编写程序求解三次样条插值函数满足 并针对下面一组具体实验数据0.250.30.390.450.530.50000.54770.62450.67080.7280求解，其中边界条件为.for i=2:n-1 A(i,i-1)=dx(i-1); A(i,i)=2*(dx(i-1)+dx(i); A(i,i+1)=dx(i); B(i,1)=3*(dy(i)/dx(i)-dy(i-1)/dx(i-1);end%端点一阶导数条件 %if flag=1 A(1,1)=2*dx(1);A(1,2)=dx(1); A(n,n-1)=dx(n-1);A(n,n)=2*dx(n-1); B(1,1)=3*(dy(1)/dx(1)-vl); B(n,1)=3*(vr-dy(n-1)/dx(n-1);end% - %端点二阶导数条件 %if flag=2 A(1,1)=2;A(n,n)=2; B(1,1)=vl;B(n,1)=vr;end% - %c=AB;for i=1:n-1 d(i)=(c(i+1)-c(i)/(3*dx(i); b(i)=dy(i)/dx(i)-dx(i)*(2*c(i)+c(i+1)/3;endclc;clear;format short g;x1=0.25 0.3 0.39 0.45 0.53;y1=0.5000 0.5477 0.6245 0.6708 0.7280;u1=0;un=0;xx1=x1(1):0.001:x1(end);yy1 b1 c1 d1=spline3(x1,y1,xx1,2,u1,un);disp(各区间上的三次样条插值系数为：)fprintf(ttb1tttc1ttttd1n);disp(b1 c1(1:end-1,1) d1);单独m文件：functionyy b c d=spline3(x,y,xx,flag,vl,vr)%样条函数为：Si(x)=yi+bi*(x-xi)+ci*(x-xi)2+di*(x-xi)3n=length(x);a=zeros(n-1,1);b=a;d=a;dx=a;dy=a;A=zeros(n);B=zeros(n,1);for i=1:n-1 a(i)=y(i); dx(i)=x(i+1)-x(i); dy(i)=y(i+1)-y(i);end mm nn=size(xx);yy=zeros(mm,nn);for i=1:mm*nn for ii=1:n-1 if xx(i)=x(ii) & xx(i)x(ii+1) j=ii; break; elseif xx(i)=x(n) j=n-1; end end yy(i)=a(j)+b(j)*(xx(i)-x(j)+c(j)*(xx(i)-x(j)2+d(j)*(xx(i)-x(j)3; endend结果如下：各区间上的三次样条插值系数为：b1c1d1 0.96966 0 -6.2652 0.92267 -0.93977 1.8813 0.79923 -0.43181 -0.46 0.74245 -0.51461 2.1442五、编写程序构造区间上的以等分结点为插值结点的Newton插值公式，假设结点数为（包括两个端点），给定相应的函数值，插值区间和等分的份数，该程序能快速计算出相应的插值公式。以，为例计算其对应的插值公式，分别取不同的值并画出原函数的图像以及插值函数的图像，观察当增大时的逼近效果.Newton插值:clc;clear;syms xa=-1;b=1;n=3;x0=a:2/(n-1):b;y0=1./(1+25.*x0.2);x1=a:0.0001:b;y1=1./(1+25.*x1.2);plot(x1,y1)hold ony=newton(x0,y0);plot(x1,subs(y,x,x1)axis(-1,1,-1,1);单独m文件：function y=newton(x0,y0)syms x Bn=length(x0);A=zeros(n,n);A(:,1)=y0; for i=2:n for j=i:n A(j,i)=(A(j,i-1)-A(j-1,i-1)/(x0(j)-x0(j-i+1); endendB(1,1)=1;for i=2:n B(i,1)=B(i-1,1)*(x-x0(i-1);endfor i=1:n C(i)=A(i,i).*B(i,1);endp=0;for i=1:n p=p+C(i);endP=expand(p);digits(6);y=vpa(P)n=3时插值结果如下：y = 1.-.961538*x2n=5时插值结果如下：y = 1.-4.27719*x2+3.31565*x4n=7时插值结果如下：y = 1.+.124928e-14*x-8.78409*x2-.284080e-14*x3+20.9574*x4-13.1349*x6n=10时插值结果如下：y= .138809e-14*x-8.26092*x2+.299257e-14*x3+30.7285*x4+.789492e-15*x7+21.6248*x8-.218329e-14*x5-44.9155*x6+.861538n=20时插值结果如下：y= .216415e-14*x+95604.8*x16-21.6235*x2+121019.*x12-58585.0*x10+.832282e-11*x9-.210106e-12*x3+327.726*x4+.390814e-11*x7+17172.5*x8+.647870e-11*x13-.921334e-12*x5-3055.32*x6+.459534e-11*x17-25671.2*x18-146791.*x14-.292343e-10*x15+.128399e-10*x11+.992681n=30时插值结果如下：y=-.577993e-8*x11+.242379e-15*x+.999614-.600443e-9*x9-.122518e7*x10-24.5956*x2-.144197e-12*x3+551.257*x4-.387656e-11*x5-9924.48*x6-.820216e-10*x7+131063.*x8+.806398e7*x12-.903419e-8*x13-.375301e8*x14+.468529e-7*x15+.123950e9*x16-.614169e-7*x17-.289894e9*x18+.339302e-7*x19+.474656e9*x20+.134317e-6*x21-.530053e9*x22-.781817e-7*x23+.383368e9*x24-.258139e-7*x25-.161437e9*x26+.411067e-8*x27+.299799e8*x28n=50时插值结果如下：y= .163155e16*x30+.288123e-14*x+.132963e9*x12-.747119e7*x10+.908504e-5*x13+.999999-.212026e-7*x9+357069.*x8+.228985e11*x16+.833425e-4*x17-.570912e16*x34-.714607*x37+.514357e-1*x43+.143133e16*x44+1.01658*x31-.19