矩阵与数值分析编程.docx
- 文档编号:23512847
- 上传时间:2023-05-17
- 格式:DOCX
- 页数:13
- 大小:67.28KB
矩阵与数值分析编程.docx
《矩阵与数值分析编程.docx》由会员分享,可在线阅读,更多相关《矩阵与数值分析编程.docx(13页珍藏版)》请在冰豆网上搜索。
矩阵与数值分析编程
矩阵与数值分析实习作业
学生班级:
01班
学生姓名:
黄晓超
任课教师:
张宏伟
所在学院:
机械学院
学生学号:
20904027
2010年1月7号
一、解线性方程组
1.分别Jacobi迭代法和Gauss-Seidel迭代法求解教材85页例题。
迭代法计算停止的条件为:
.
解:
求线性方程组,
3.4336-0.523800.67105-0.15270x1-1.0
-0.523803.28326-03.73051-0.26890x21.5
0.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
%Detailedexplanationgoeshere
clc;
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);
while1
fori=1:
n
y(i)=b(i);
forj=1:
n
ifj~=i
y(i)=y(i)-A(i,j)*x(j);
end
end
y(i)=y(i)/A(i,i);
end
ifnorm(y-x,inf) break; end x=y; k=k+1; end y Jacobi法程序运算结果: y= -0.3941 0.5058 0.7612 -0.7030 2、Gauss-Seidel迭代法编译程序: function[output_args]=Untitled2(input_args)%G-S迭代法 %UNTITLED2Summaryofthisfunctiongoeshere %Detailedexplanationgoeshere clc; 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); while1 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(djwcX break; end end X Gauss-Seidel迭代法运行的解 X= -0.3941 0.5058 0.7612 -0.7030 3、分析: 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.5 clc; 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; end x=Xk1%x的值即为所求 i Newton迭代法运行结果 x= 0.5671 i= 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; end x=Xk3 i 割线法 运行结果 x= 0.5671 i= 5 3、比较分析: 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))>=e Q=0; a=0; whilea<=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; forj=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; end H1(i+1)%存放算法求得的解 function[y]=f(x) y=2/sqrt(pi)*exp(-x^2); Romberg方法的程序计算结果 ans= 0.8427 四、插值与逼近 2.已知函数值 0 1 2 3 4 5 6 7 8 9 10 0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29 和边界条件: . 求三次样条插值函数 并画出其图形. 解: 编译程序: : 在工作空间输入: clear closeall x=0: 10; y=[00.791.532.192.713.033.272.893.063.193.29]; %求解mi程序 gk=[]; fork=1: 9 gk(k)=3*(y(k+2)-y(k))/2; end m0=0.8; m10=0.2; A=[20.50000000 0.520.5000000 00.520.500000 000.520.50000 0000.520.5000 00000.520.500 000000.520.50 0000000.520.5 00000000.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),gridon title('三次样条插值图(第一类边界条件)') xlabel('x') ylabel('y') 结果: 数据结果: m= 0.77148596314996 0.70405614740014 0.61228944724946 0.38678606360200 0.36056629834254 -0.14905125697216 -0.18436127045388 0.25649633878770 0.05837591530307 图形如下图: 五、常微分方程数值解法 用Euler法与改进的Euler法求解 取步长 计算,画出数值解的图形并与精确解 相比较。 解: 编译程序 function[output_args]=Untitled1(input_args) %UNTITLED1Summaryofthisfunctiongoeshere %Detailedexplanationgoeshere clc; 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: N Y(j+1)=Y(j)+h*(Y(j)-T(j)*Y(j)^2); end plot(T,Y,'g');gridon; holdon; %Euler改进法 forj=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); end plot(T,Y,'r'); holdon; %精确解 Y=1./(T-1+2*exp(-T)); plot(T,Y); xlabel('x'); ylabel('f(x)'); title('Euler(绿色)、Euler改进法(红色)及精确解(蓝色)的比较'); 分析: Euler法收敛阶是1阶,改进的Euler算法精度比较高更加逼近真实值。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 矩阵 数值 分析 编程
![提示](https://static.bdocx.com/images/bang_tan.gif)