矩阵与数值分析编程

矩阵与数值分析实习作业学生班级：01班学生姓名：黄晓超任课教师：张宏伟所在学院：机械学院学生学号：209040272010年1月7号一、解线性方程组1．分别Jacobi迭代法和Gauss-Seidel迭代法求解教材85页例题。迭代法计算停止的条件为：．解：求线性方程组,3.4336-0.523800.67105-0.15270x1-1.0-0.523803.28326-03.73051-0.26890x21.50.67105-0.730514.02612-0.009835x32.5-0.15272-0.268900.018352.75702x4-2.0使用MATLAB编译程序1、Jacobi迭代法程序：function[output_args]=Untitled1(input_args)%Jacobi法%UNTITLED1Summaryofthisfunctiongoeshere%Detailedexplanationgoeshereclc;clear;A=[3.4336-0.523800.67105-0.15270;-0.523803.28326-0.73051-0.26890;0.67105-0.730514.02612-0.09835;-0.15272-0.268900.018352.75702];b=[-1.01.52.5-2.0];delta=10^(-6);%误差n=length(A);k=0;x=zeros(n,1);y=zeros(n,1);while1fori=1:ny(i)=b(i);forj=1:nifj~=iy(i)=y(i)-A(i,j)*x(j);endendy(i)=y(i)/A(i,i);endifnorm(y-x,inf)<deltabreak;endx=y;k=k+1;endyJacobi法程序运算结果：y=-0.39410.50580.7612-0.70302、Gauss-Seidel迭代法编译程序：function[output_args]=Untitled2(input_args)%G-S迭代法%UNTITLED2Summaryofthisfunctiongoeshere%Detailedexplanationgoeshereclc;clear;A=[3.4336-0.523800.67105-0.15270;-0.523803.28326-0.73051-0.26890;0.67105-0.730514.02612-0.09835;-0.15272-0.268900.018352.75702];b=[-1.01.52.5-2.0];delta=10^(-6);%误差X=zeros(4,1);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);jX=A\b';[nm]=size(A);iD=inv(D-L);B2=iD*U;H=eig(B2);mH=norm(H,inf);while1iD=inv(D-L);B2=iD*U;f2=iD*b';X1=B2*X+f2;X=X1;djwcX=norm(X1-jX,inf);xdwcX=djwcX/(norm(X,inf)+eps);if(djwcX<delta)|(xdwcX<delta)break;endendXGauss-Seidel迭代法运行的解X=-0.39410.50580.7612-0.70303、分析：Jacobi迭代法称简单迭代法，程序编译比较简单，但是在迭代过程中，对已经算出来的信息未加充分利用，但是该方法算出来的值比较精确，Gauss-Seidel迭代法占用内存小，收敛快，但是精确度稍差。二、非线性方程的迭代解法1．用Newton迭代法、割线法求方程在附近的正．计算停止的条件为：；解：首先使用使用Newton迭代法1、MATLAB的编译程序如下：function[output_args]=Untitled1(input_args)%UNTITLED1Summaryofthisfunctiongoeshere%Detailedexplanationgoeshere%f(x)=x*exp(x)-1%f'(x)=exp(x)+x*exp(x)%φ'(x)=x-f(x)/f'(x)%迭代初值为x=0.5clc;clear;i=1;Xk=0.5;Xk1=Xk-(Xk*exp(Xk)-1)/(exp(Xk)+Xk*exp(Xk));whileabs(Xk1-Xk)>=10^(-6)Xk=Xk1;Xk1=Xk-(Xk*exp(Xk)-1)/(exp(Xk)+Xk*exp(Xk));i=i+1;endx=Xk1%x的值即为所求iNewton迭代法运行结果x=0.5671i=4再运用割线法割线法2、MATLAB的编译程序如下function[output_args]=Untitled2(input_args)%UNTITLED2Summaryofthisfunctiongoeshere%Detailedexplanationgoeshere%迭代初值x1=0.5,x2=0.4%迭代公式为Xk+1=Xk-f(Xk)/((f(Xk)-f(Xk-1))/(Xk-Xk-1))clc;clear;i=1;Xk1=0.5;Xk2=0.4;Xk3=Xk2-(Xk2*exp(Xk2)-1)*(Xk2-Xk1)/((Xk2*exp(Xk2)-1)-(Xk1*exp(Xk1)-1));whileabs(Xk3-Xk2)>=10^(-6)Xk1=Xk2;Xk2=Xk3;Xk3=Xk2-(Xk2*exp(Xk2)-1)*(Xk2-Xk1)/((Xk2*exp(Xk2)-1)-(Xk1*exp(Xk1)-1));i=i+1;endx=Xk3i割线法运行结果x=0.5671i=53、比较分析：Newton迭代法是二阶收敛，割线法是在前者基础上变化来的，收敛阶1.618，在该题中前者的迭代步数为4，后者为5，可见Newton迭代法收敛比较快。三、数值积分用Romberg方法求的近似值，要求误差不超过.解：Romberg的编译程序：function[output_args]=Untitled1(input_args)%UNTITLED1Summaryofthisfunctiongoeshere%Detailedexplanationgoeshere%求解2/sqre(pi)*|exp(-x^2)(01)的值。clc;clear;e=10^(-5);%误差i=1;%循环控制变量H0(i)=(f(0)+f(1))/2;H1(i)=H0(i)/2+f(0.5)/2;H1(i+1)=H1(i)*4/3-H0(i)/3;%两层循环whileabs(H1(i+1)-H0(i))>=eQ=0;a=0;whilea<=1a=a+1/2^(i+1);Q=Q+f(a);a=a+1/2^(i+1);endQ=Q/2^(i+1);H2(1)=H1(1)/2+Q;forj=2:i+2H2(j)=4^(j-1)/(4^(j-1)-1)*H2(j-1)-1/(4^(j-1)-1)*H1(j-1);endH0=H1;H1=H2;i=i+1;endH1(i+1)%存放算法求得的解function[y]=f(x)y=2/sqrt(pi)*exp(-x^2);Romberg方法的程序计算结果ans=0.8427四、插值与逼近2．已知函数值01234567891000.791.532.192.713.033.272.893.063.193.29和边界条件：．求三次样条插值函数并画出其图形．解：编译程序:：在工作空间输入：clearcloseallx=0:10;y=[00.791.532.192.713.033.272.893.063.193.29];%求解mi程序gk=[];fork=1:9gk(k)=3*(y(k+2)-y(k))/2;endm0=0.8;m10=0.2;A=[20.500000000.520.500000000.520.500000000.520.500000000.520.500000000.520.500000000.520.500000000.520.500000000.52];b1=gk(1)-0.5*m0;b9=gk(9)-0.5*m10;b=[b1gk(2)gk(3)gk(4)gk(5)gk(6)gk(7)gk(8)b9]';m=A\b%图形程序pp=csape(x,y,'complete',[0.8,0.2]);%第一类边界条件调用函数xi=0:0.25:10;yi=ppval(pp,xi);plot(x,y,'o',xi,yi),gridontitle('三次样条插值图（第一类边界条件）')xlabel('x')ylabel('y')结果：数据结果：m=0.771485963149960.704056147400140.612289447249460.386786063602000.36056629834254-0.14905125697216-0.184361270453880.256496338787700.05837591530307图形如下图：五、常微分方程数值解法用Euler法与改进的Euler法求解取步长计算，画出数值解的图形并与精确解相比较。解：编译程序function[output_args]=Untitled1(input_args)%UNTITLED1Summaryofthisfunctiongoeshere%Detailedexplanationgoeshereclc;clear;h=0.1;a=0;b=1;ya=1;N=(b-a)/h;T=linspace(a,b,N+1);Y=zeros(1,N+1);Y(1)=ya;%Euler法forj=1:NY(j+1)=Y(j)+h*(Y(j)-T(j)*Y(j)^2);endplot(T,Y,'g');gridon;holdon;%Euler改进法forj=1:NY(j+1)=Y(j)+(h/2)