级大连理工矩阵实验报告.docx
- 文档编号:18225166
- 上传时间:2023-08-14
- 格式:DOCX
- 页数:32
- 大小:101.45KB
级大连理工矩阵实验报告.docx
《级大连理工矩阵实验报告.docx》由会员分享,可在线阅读,更多相关《级大连理工矩阵实验报告.docx(32页珍藏版)》请在冰点文库上搜索。
级大连理工矩阵实验报告
2009级《矩阵与数值分析》实验报告
2010年1月5号
2009级工科硕士研究生
《矩阵与数值分析》课程数值实验题目
一、解线性方程组
1.分别Jacobi迭代法和Gauss-Seidel迭代法求解教材85页例题。
迭代法计算停止的条件为:
.
解:
首先用Jacobi迭代法进行编程,具体程序如下:
#include"stdafx.h"
#include
#include
usingnamespacestd;
voidmain()
{
doublea[4][4]={{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[4]={-1.0,1.5,2.5,-2.0},x[4]={0},x1[4],r=1.0e-6,t,s;
inti,j,c,d,e=0;
cout.precision(8);
do
{
for(c=0;c<4;c++)
{
x1[c]=x[c];
}
for(i=0;i<4;i++)
{
s=0.0;
for(j=0;j<4;j++)
{
if(j!
=i)
S+=a[i][j]*x[j];
}
x[i]=(b[i]-s)/a[i][i];
}
t=fabs(x[0]-x1[0]);
for(d=0;d<4;d++)
{
cout< if(t t=fabs(x[d]-x1[d]); } cout< e++; if(t break; }while (1); cout<<"迭代次数为: "< } 运行结果为: -0.291239520.410399650.74395133-0.7064777 -0.405446280.499844440.76195792-0.70420001 -0.395219210.505668960.76136578-0.70306148 -0.394164320.505798750.76124132-0.70298956 -0.394116990.50578450.76123261-0.70298827 -0.394117410.50578260.76123236-0.70298847 -0.394117660.505782490.76123238-0.7029885 迭代次数为: 7 再用Gauss-Seidel迭代法进行编程,具体程序如下: #include"stdafx.h" #include #include usingnamespacestd; voidmain() { doublea[4][4]={{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[4]={-1.0,1.5,2.5,-2.0},x[4]={0},x1[4],r=1.0e-6,t,s; inti,j,c,d,e=0; cout.precision(8); do { for(c=0;c<4;c++) { x1[c]=x[c]; } for(i=0;i<4;i++) { s=0.0; for(j=0;j<4;j++) { if(j! =i) S+=a[i][j]*x1[j]; } x[i]=(b[i]-s)/a[i][i]; } t=fabs(x[0]-x1[0]); for(d=0;d<4;d++) { cout< if(t t=fabs(x[d]-x1[d]); } cout< e++; if(t break; }while (1); cout<<"迭代次数为: "< } 运行结果为: -0.291239520.4568630.62094523-0.72542093 -0.375160830.489144860.73466119-0.7011273 -0.391380020.503047270.75509946-0.70338429 -0.393353950.504822280.76027013-0.70306281 -0.394079410.505684150.76092905-0.70303344 -0.39407540.505717430.76120706-0.70299395 -0.39412290.505783160.7612134-0.70299234 -0.394114040.505777120.76123328-0.7029886 -0.394118680.505783260.7612308-0.70298883 -0.394117270.505781950.76123268-0.70298847 -0.394117820.505782630.76123222-0.70298853 迭代次数为: 11 分析: 结果为x=[-0.394117820.505782630.76123222-0.70298853]';Jacobi迭代法迭代次数为11,而Gauss-Seidel迭代法次数为7,迭代的次数更少。 说明Gauss-Seidel迭代法收敛性更好。 2.用Gauss列主元消去法、QR方法求解如下方程组: 解: 首先用Gauss列主元消去法来编程求解,具体程序如下: %Gauss列主元消去法 #include"stdafx.h" #include usingnamespacestd; voidmain() { floata[4][4]={{2,4,3,1},{8,2,0,0},{5,0,4,0},{9,0,0,5}}; floatb[4]={12,6,23,16}; intt,s,i,j,k; floatx[4]; for(i=0;i<4;i++) { t=a[i][i]; for(j=i;j<4;j++) { if(a[j][i]>t) { s=j; t=a[j][i]; } } floatc; for(k=i;k<4;k++) { c=a[i][k]; a[i][k]=a[s][k]; a[s][k]=c; } c=b[i]; b[i]=b[s]; b[s]=c; for(j=i+1;j<4;j++) { c=a[j][i]/a[i][i]; for(k=i;k<4;k++) { a[j][k]-=c*a[i][k]; } b[j]-=c*b[i]; } } for(i=3;i>=0;i--) { for(j=i+1;j<4;j++) { b[i]-=x[j]*a[i][j]; } x[i]=b[i]/a[i][i]; } cout<<"Theansweris: "< for(i=0;i<4;i++) { cout<<"x"< } cout< } 运行结果: Theansweris: x1=1.04604x2=-1.18414x3=4.44246x4=1.31714 再用QR方法来求解,程序如下: A=[2,4,3,1;8,2,0,0;5,0,4,0;9,0,0,5]; b=[12,6,23,16]'; [Q,R]=qr(A) x=R\(Q\b) 运行结果: Q=-0.15160.9116-0.20050.3253 -0.60650.21950.4010-0.6505 -0.3790-0.1688-0.8765-0.2439 -0.6823-0.30390.17500.6415 R=-13.1909-1.8194-1.9711-3.5631 04.08532.0595-0.6077 00-4.10770.6747 0003.5327 x=1.0460 -1.1841 4.4425 1.3171 二、非线性方程的迭代解法 1.用Newton迭代法、割线法求方程 在 附近的正.计算停止的条件为: ; 解: 用Newton迭代法来求解,程序如下: x0=0.5; i=0; deta=0.1; whiledeta>=1.0e-6 f=x0*exp(x0)-1; f1=exp(x0)+x0*exp(x0); x=x0-f/f1 deta=abs(x-x0) x0=x; i=i+1; end '迭代次数' i '迭代结果' x 运行结果: x=0.5710 deta=0.0710 x=0.5672 deta=0.0039 x=0.5671 deta=1.2278e-005 x=0.5671 deta=1.2348e-010 ans= 迭代次数 i=4 ans= 迭代结果 x=0.5671 用Newton割线法来求解,再取一个初值设x1=1.4,程序如下: x0=0.5; x1=0.4; i=0; deta=0.1; whiledeta>=1.0e-6 f0=x0*exp(x0)-1; f1=x1*exp(x1)-1; x=x1-f1/(f1-f0)*(x1-x0) deta=abs(x-x0) x0=x1; x1=x; i=i+1; end '迭代次数' i '迭代结果' x 运行结果如下: x=0.5772 deta=0.0772 x=0.5657 deta=0.1657 x=0.5671 deta=0.0100 x=0.5671 deta=0.0014 x=0.5671 deta=1.1854e-005 x=0.5671 deta=1.4080e-008 ans= 迭代次数 i=6 ans= 迭代结果 x=0.5671 分析: 运行结果得到x=0.5671,Newton迭代法迭代次数为4,Newton割线法迭代次数为6,说明本题用Newton迭代法效果更好。 2.利用迭代法求多项式 的所有零点。 解: 利用Newton迭代法求解每个区间上的根,程序如下: forx=-100: 0.1: 100 f=x^5-13.7*x^4+14.66*x^3+213.092*x^2-154.2072*x-631.9296; x=x+0.01; f1=x^5-13.7*x^4+14.66*x^3+213.092*x^2-154.2072*x-631.9296; iff*f1<0 x0=x; x1=x+0.1; i=0; deta=0.1; whiledeta>=1.0e-8 f0=x0^5-13.7*x0^4+14.66*x0^3+213.092*x0^2-154.2072*x0-631.9296; f1=x1^5-13.7*x1^4+14.66*x1^3+213.092*x1^2-154.2072*x1-631.9296; x=x1-f1/(f1-f0)*(x1-x0); deta=abs(x-x0); x0=x1; x1=x; i=i+1; end '迭代次数' i '迭代结果' x end end 运行结果如下: ans= 迭代次数 i=1 ans= 迭代结果 x=-3 ans= 迭代次数 i=6 ans= 迭代结果 x=-1.6000 ans= 迭代次数 i=6 ans= 迭代结果 x=2.3000 ans= 迭代次数 i=7 ans= 迭代结果 x=5.4000 ans= 迭代次数 i=7 ans= 迭代结果 x=10.6000 分析: 首先判断函数的大致的根的范围,再用迭代法莱求解,最后求出根分别为-3,-1.6,2.3,5.4,10.6. 三、数值积分 用Romberg方法求 的近似值,要求误差不超过 . 解: 根据要求进行编程,具体程序如下: %RombergIntegration clc; a=0;b=1; h=b-a; symsx; error=1; i=2; r=[]; f=2./sqrt(pi)*exp(-x.^2); r(1,1)=h*(subs(f,a)+subs(f,b))/2 whileerror>1.0e-5 sum=0; fork=1: 2^(i-2) sum=sum+subs(f,a+(k-0.5)*h); end; r(2,1)=(r(1,1)+h*sum)/2; forj=2: i r(2,j)=r(2,j-1)+(r(2,j-1)-r(1,j-1))/(4^(j-1)-1); end; r(2,: ) error=abs(r(2,i)-r(1,i-1)); h=h/2; i=i+1; r(1,: )=r(2,: ); end 运行结果: r=0.77174 ans=0.825260.8431 ans=0.838370.842740.84271 ans=0.841620.84270.84270.8427 ans=0.842430.84270.84270.84270.8427 分析: 最后得到在误差范围内的根为0.8427。 四、插值与逼近 1.给定 上的函数 ,请做如下的插值逼近: 构造等距节点分别取 , , 的Lagrange插值多项式; 解: 根据要求进行编程,多项式的变量用t来表示的,程序如下: %拉格朗日插值等分5个节点 X1=linspace(-1,1,5); fori=1: 5 Y1(i)=1./(1+25*X1(i).^2); end p1=0; symst; fori=1: 5 L=1; forj=1: 5 ifj~=i L=L*(t-X1(j))/(X1(i)-X1(j)); end end p1=p1+L*Y1(i); end p1=simplify(p1); p1=vpa(p1,4) %拉格朗日插值等分8个节点 X2=linspace(-1,1,8); fori=1: 8 Y2(i)=1./(1+25*X2(i).^2); end p2=0; symst; fori=1: 8 L=1; forj=1: 8 ifj~=i L=L*(t-X2(j))/(X2(i)-X2(j)); end end p2=p2+L*Y2(i); end p2=simplify(p2); p2=vpa(p2,4) %拉格朗日插值等分10个节点 X3=linspace(-1,1,10); fori=1: 10 Y3(i)=1./(1+25*X3(i).^2); end p3=0; symst; fori=1: 10 L=1; forj=1: 10 ifj~=i L=L*(t-X3(j))/(X3(i)-X3(j)); end end p3=p3+L*Y3(i); end p3=simplify(p3); p3=vpa(p3,4) %曲线图 figure;holdon; t=-1: 0.001: 1; plot(t,subs(p1,t),'b');gridon; plot(t,subs(p2,t),'r'); plot(t,subs(p3,t),'g'); title('函数曲线'); xlabel('自变量');ylabel('函数值'); 运行结果如下: p1=3.978*t^4-4.939*t^2+1. p2=-4.617*t^2+9.077*t^4-5.174*t^6+.7527 p3=-8.263*t^2+30.74*t^4-44.93*t^6+21.63*t^8+.8616 分析: 利用Lagrange插值求多项式时根据等分的节点数不同得到的曲线的拟合的误差也不同,在本题中,节点增大,拟合的越接近,误差越小。 (2)构造分段线性取 的Lagrange插值多项式; 解: 根据分段线性插值公式进行编程,程序如下: %拉格朗日插值等分10个节点,分段线性插值 clc n=10; X=linspace(-1,1,n); fori=1: n Y(i)=1./(1+25*X(i).^2); end symsx; L=0; fori=1: 9 L=Y(i)*(x-X(i+1))/(X(i)-X(i+1))+Y(i+1)*(x-X(i))/(X(i+1)-X(i)); L=vpa(L,4) End 运行结果得到各个区间上的插值多项式: L=.1060*x+.1445 L=.2372*x+.2465 L=.6749*x+.4897 L=2.248*x+1.014 L=.7642 L=-2.248*x+1.014 L=-.6749*x+.4897 L=-.2372*x+.2465 L=-.1060*x+.1445 运行程序: %分段线性插值,等分10个节点,生成原函数的曲线和插值后的曲线 clc n=10; x0=linspace(-1,1,n); fori=1: n y0(i)=1./(1+25*x0(i).^2); end forx=-1: 0.001: 1 forj=1: length(x) fori=1: n-1 if(x>=x0(i))&&(x<=x0(i+1)) y=(x-x0(i+1))/(x0(i)-x0(i+1))*y0(i)+(x-x0(i))/(x0(i+1)-x0(i))*y0(i+1); else continue; end end end holdon; plot(x,y,'r'); plot(x,1./(1+25*x.^2)); end 运行结果: 分析: 利用分段线性Lagrange插值来拟合曲线时可以很好的接近原曲线,拟合的误差很小。 比节点法得到的结果要好。 取Chebyshev多项式 的零点: , 作插值节点构造 的插值多项式 和上述的插值多项式均要求画出曲线图形(用不同的线型或颜色表示不同的曲线)。 解: 根据要求进行编程得到程序如下: %利用Chebyshev多项式的零点作为插值节点 clc; n=10; fork=1: n X(k)=cos((2*k-1)/(2*n)*pi); Y(k)=1./(1+25*X(k).^2); end p=0; symst; L=1; fori=1: n L=1; forj=1: n ifj~=i L=L*(t-X(j))/(X(i)-X(j)); end end p=p+L*Y(i); end p=simplify(p); p=vpa(p,4) t=-1: 0.001: 1; figure;holdon; plot(t,subs(1./(1+25*t.^2),t),'black');gridon; plot(t,subs(p,t)); title('函数曲线'); xlabel('自变量');ylabel('函数值'); 运行结果: p= .5612e-15*t+.4767e-14*t^9+.7308-4.812*t^2-.4795e-14*t^3+12.62*t^4+.1288e-13*t^5-14.00*t^6-.1345e-13*t^7+5.513*t^8 运行曲线图: 分析: 利用Chebyshev多项式的零点来作插值节点比利用等分区间上的节点作插值节点误差要小,拟合的曲线更接近真实的曲线。 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 和边界条件: . 求三次样条插值函数 并画出其图形. 解: 根据要求得到程序如下: %三次样条插值 closeall clearall clc X=[0,1,2,3,4,5,6,7,8,9,10]; Y=[0,0.79,1.53,2.19,2.71,3.03,3.27,2.89,3.06,3.19,3.29]; m (1)=0.8;m(11)=0.2;h=1;lmt=0.5;mu=0.5;n=11;g=[]; fork=2: 10 g(k)=1.5*(Y(k+1)-Y(k-1)); end M=zeros(9,1);g (2)=g (2)-0.4;g(10)=g(10)-0.1;A=zeros(9,9); fork=1: 9 A(k,k)=2; G(k)=g(k+1); end fork=1: 8 A(k,k+1)=mu; end fork=2: 9 A(k,k-1)=lmt; end M=inv(A)*G';m(2: 10)=M(1: 9); symsx; ans=[]; fork=1: 10 s=(1+2*(x-X(k)))*(x-X(k+1))^2*Y(k)+(1-2*(x-X(k+1)))*(x-X(k))^2*Y(k+1)+(x-X(
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 大连理工 矩阵 实验 报告
![提示](https://static.bingdoc.com/images/bang_tan.gif)