矩阵与数值分析编程.docx

1、矩阵与数值分析编程矩阵与数值分析实习作业学生班级：01班学生姓名：黄晓超任课教师：张宏伟所在学院：机械学院学生学号：209040272010年1月7号一、解线性方程组1分别Jacobi迭代法和Gauss-Seidel迭代法求解教材85页例题。迭代法计算停止的条件为：解：求线性方程组,3.4336 -0.52380 0.67105 -0.15270 x1 -1.0-0.52380 3.28326 -03.73051 -0.26890 x2 1.50.67105 -0.73051 4.02612 -0.009835 x3 2.5-0.15272 -0.26890 0.01835 2.75702 x

2、4 -2.0使用MATLAB编译程序1、Jacobi迭代法程序：function output_args = Untitled1( input_args )%Jacobi法%UNTITLED1 Summary of this function goes here% Detailed explanation goes hereclc;clear;A=3.4336 -0.52380 0.67105 -0.15270;-0.52380 3.28326 -0.73051 -0.26890; 0.67105 -0.73051 4.02612 -0.09835;-0.15272 -0.26890 0.01

3、835 2.75702;b=-1.0 1.5 2.5 -2.0;delta=10(-6);%误差n=length(A);k=0;x=zeros(n,1);y=zeros(n,1);while 1 for i=1:n y(i)=b(i); for j=1:n if j=i y(i)=y(i)-A(i,j)*x(j); end end y(i)=y(i)/A(i,i); end if norm(y-x,inf)delta break; end x=y; k=k+1;endyJacobi法程序运算结果：y = -0.3941 0.5058 0.7612 -0.70302、Gauss-Seidel迭代

4、法编译程序：function output_args = Untitled2( input_args )%G-S迭代法%UNTITLED2 Summary of this function goes here% Detailed explanation goes hereclc;clear;A=3.4336 -0.52380 0.67105 -0.15270;-0.52380 3.28326 -0.73051 -0.26890; 0.67105 -0.73051 4.02612 -0.09835;-0.15272 -0.26890 0.01835 2.75702;b=-1.0 1.5 2.5

5、-2.0;delta=10(-6);%误差X=zeros(4,1);D=diag(diag(A);U=-triu(A,1);L=-tril(A,-1); jX=Ab;n m=size(A); iD=inv(D-L); B2=iD*U; H=eig(B2);mH=norm(H,inf);while 1 iD=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 (djwcXdelta)|(xdwcX=10(-6) Xk=Xk1; Xk1=Xk

6、-(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 )%UNTITLED2 Summary of this function goes here% Detailed explanation goes here%迭代初值x1=0.5,x2=0.4%迭代公式为Xk+1=Xk-f(Xk)/(f(Xk)-f(Xk-1)/(Xk-Xk-1)cl

7、c;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);while abs(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，在该题中前者

8、的迭代步数为4，后者为5，可见Newton迭代法收敛比较快。三、数值积分用Romberg方法求的近似值，要求误差不超过.解：Romberg的编译程序：function output_args = Untitled1( input_args )%UNTITLED1 Summary of this function goes here% Detailed explanation goes here%求解2/sqre(pi)*|exp(-x2) (0 1)的值。clc;clear;e=10(-5);%误差i=1;%循环控制变量H0(i)=(f(0)+f(1)/2;H1(i)=H0(i)/2+f(0.

9、5)/2;H1(i+1)=H1(i)*4/3-H0(i)/3;%两层循环while abs(H1(i+1)-H0(i)=e Q=0; a=0; while a=1 a=a+1/2(i+1); Q=Q+f(a); a=a+1/2(i+1); end Q=Q/2(i+1); H2(1)=H1(1)/2+Q; for j=2:i+2 H2(j)=4(j-1)/(4(j-1)-1)*H2(j-1)-1/(4(j-1)-1)*H1(j-1); end H0=H1; H1=H2; i=i+1;endH1(i+1)%存放算法求得的解function y = f(x)y=2/sqrt(pi)*exp(-x2)

10、;Romberg方法的程序计算结果ans =0.8427四、插值与逼近2已知函数值01234567891000.791.532.192.713.033.272.893.063.193.29和边界条件：求三次样条插值函数并画出其图形解 ：编译程序: ：在工作空间输入： clear close all x=0:10; y=0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29; %求解mi程序 gk=; for k=1:9 gk(k)=3*(y(k+2)-y(k)/2; end m0=0.8; m10=0.2; A= 2 0.5 0 0 0 0 0

11、0 0 0.5 2 0.5 0 0 0 0 0 0 0 0.5 2 0.5 0 0 0 0 0 0 0 0.5 2 0.5 0 0 0 0 0 0 0 0.5 2 0.5 0 0 0 0 0 0 0 0.5 2 0.5 0 0 0 0 0 0 0 0.5 2 0.5 0 0 0 0 0 0 0 0.5 2 0.5 0 0 0 0 0 0 0 0.5 2 ; b1=gk(1)-0.5*m0; b9=gk(9)-0.5*m10; b=b1 gk(2) gk(3) gk(4) gk(5) gk(6) gk(7) gk(8) b9; m=Ab %图形程序 pp=csape(x,y,complete,0

12、.8,0.2); %第一类边界条件调用函数 xi=0:0.25:10; yi=ppval(pp,xi); plot(x,y,o,xi,yi),grid on title(三次样条插值图（第一类边界条件）) xlabel(x) ylabel(y) 结果： 数据结果： m = 0.77148596314996 0.70405614740014 0.61228944724946 0.38678606360200 0.36056629834254 -0.14905125697216 -0.18436127045388 0.25649633878770 0.05837591530307 图形如下图：五、

13、常微分方程数值解法用Euler法与改进的Euler法求解取步长计算，画出数值解的图形并与精确解相比较。解：编译程序function output_args = Untitled1( input_args )%UNTITLED1 Summary of this function goes here% Detailed explanation goes hereclc;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法for j=1:N Y(j+1)=Y(j)+h*(Y(j)-T

14、(j)*Y(j)2);endplot(T,Y,g);grid on;hold on;%Euler改进法for j=1:N Y(j+1)=Y(j)+(h/2)*(Y(j)-T(j)*Y(j)2+ (Y(j)+h*(Y(j)-T(j)*Y(j)2)-T(j+1)*(Y(j)+h*(Y(j)-T(j)*Y(j)2)2 );endplot(T,Y,r);hold on;%精确解Y=1./(T-1+2*exp(-T);plot(T,Y);xlabel(x);ylabel(f(x);title(Euler(绿色)、Euler改进法(红色)及精确解(蓝色)的比较);分析：Euler法收敛阶是1阶，改进的Euler算法精度比较高更加逼近真实值。