矩陣與數值分析編程

矩陣與數值分析實習作業(yè)學生班級：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;%循環(huán)控制變量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;%兩層循環(huán)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)