矩阵报告.docx
- 文档编号:29995720
- 上传时间:2023-08-04
- 格式:DOCX
- 页数:17
- 大小:164.43KB
矩阵报告.docx
《矩阵报告.docx》由会员分享,可在线阅读,更多相关《矩阵报告.docx(17页珍藏版)》请在冰豆网上搜索。
矩阵报告
矩阵与数值分析实习报告
00xx00xx
00x
00xxx
一、解线性方程组
对于线性方程组
(1)用直接法求解;
(2)用Jacobi迭代法求解;
(3)分别取
,用SOR方法求解.比较迭代结果(与精确解比较).
解:
(1)直接法求解
程序:
(direct.m)
A=[-4,1,1,1;1,-4,1,1;1,1,-4,1;1,1,1,-4];%系数矩阵
b=diag(ones(4));
x=A\b%直接法求解并输出
运行结果:
(2)Jacobi迭代法求解:
迭代公式:
程序:
(jacobi.m)
clearall;
formatlongg;
A=[-4,1,1,1;1,-4,1,1;1,1,-4,1;1,1,1,-4];%系数矩阵
b=[1111];
[M,N]=size(A);
x=zeros(M,1);%给方程设置初值
ep=1e-6;%设定迭代终止条件
D=diag(diag(A));
L=tril(tril(A,-1),-1);
U=triu(triu(A,1),1);
B=-inv(D)*(L+U);%求得Jacobi法迭代矩阵B
d=inv(D)*b;%求d
lamada=eig(B);
lamada=abs(lamada);%求谱半径
lamada=max(lamada);
if(lamada>1)%判断迭代方法是否收敛,决定是否运行
error('方程发散');
else
Xk=B*x+d;
n=1;%n迭代次数
while1
e=Xk-x;
e=abs(e);
lamada=max(e);
iflamada break; else%根据收敛条件决定迭代次数 x=Xk; Xk=B*x+d; n=n+1; end end end fprintf('方程组的解x=');%输出方程组的解 x fprintf('迭代次数n=');%输出迭代次数 n 运行结果: (3)SOR方法求解 基本公式: 程序: (SOR.m) formatlongg; A=[-4111;1-411;11-41;111-4];%系数矩阵 b=[1111]'; x0=[0000]';%迭代初始值 w=0.75;%松弛因子 D=diag(diag(A));%下三角矩阵 U=-triu(A,1);%上三角矩阵 L=-tril(A,-1); lw=(D-w*L)\((1-w)*D+w*U); f=((D-w*L)\b)*w; y=lw*x0+f;%迭代公式 n=1; whilenorm(y-x0)>=1.0e-6&n<=1000%选取1.0e-6作为迭代精度 x0=y; y=lw*x0+f; n=n+1; end fprintf('方程组的解y=');%输出方程组的解 y fprintf('迭代次数n=');%输出迭代次数 n 运行结果(ω=0.75的仿真结果): 分别将程序中的松弛因子ω的值修改为0.75,1.0,1.25,1.5得出结果。 迭代结果比较如下表: X1 X2 X3 X4 精确解 -1 -1 -1 -1 SOR法 解 ω=0.75 -0.999998701967 -0.999999377426 -0.99999990765 -0.999999990151 ω=1.0 -0.999998799216 -0.999999459060 -1.00000003654 -0.999999831038 ω=1.25 -0.999998889180 -0.999999529990 -1.00000002055 -1.000000176043 ω=1.5 -0.999998972403 -0.999999591619 -0.99999996115 -0.999999895439 二、三次样条插值 已知函数值 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 和边界条件: . 求三次样条插值函数 并画出其图形. 解: 程序: (chazhi.m) formatlongg; xk=[012345678910];%给定节点 y=[00.791.532.192.713.033.272.893.063.193.29];%节点处的值 h=zeros(1,10);%分别设置初始状态为0向量 c=zeros(1,9);d=zeros(1,9);g=zeros(1,9);b=zeros(1,9);%设置初始状态均为0向量 m2=zeros(9,1);m1=zeros(9,1); m0=0.8;m10=0.2;%已知的两个边界条件 fori=1: 1: 10 h(i)=xk(i+1)-xk(i);%算出hk end fori=1: 1: 9 c(i)=h(i+1)/(h(i+1)+h(i)); d(i)=h(i)/(h(i+1)+h(i)); end fori=1: 1: 9 g(i)=3*(d(i)*(y(i+2)-y(i+1))/h(i)+c(i)*(y(i+1)-y(i))/h(i)); end v=[222222222];%对角元均为2 A=diag(v);%生成矩阵 fori=1: 1: 8 A(i,i+1)=d(i); A(i+1,i)=c(i+1); b(i)=g(i); end b (1)=g (1)-c (1)*m0; b(9)=g(9)-d(9)*m10; [L,U]=lu(A);%将A进行LU分解,解方程,得出m向量 m1=L\b';%LU法求解 m2=U\m1; m (1)=m0;m(11)=m10; fori=2: 1: 10 m(i)=m2(i-1); end fori=1: 1: 10;%计算S(x)的格式 x=i-1: 0.1: i;sx(i,: )=(1-2.*(x-xk(i))./(-h(i))).*(((x-xk(i+1))./(-h(i))).^ (2)).*y(i)+... (1-2.*(x-xk(i+1))./h(i)).*(((x-xk(i))./h(i)).^ (2)).*y(i+1)+... (x-xk(i)).*(((x-xk(i+1))./-h(i)).^ (2)).*m(i)+... (x-xk(i+1)).*(((x-xk(i))./h(i)).^ (2)).*m(i+1); end sx1=zeros(1,101); x1=0: 0.1: 10; forj=1: 10; fori=1: 11; sx1(1,i+10*(j-1))=sx(j,i); end end plot(x1,sx1);%画出图形 sx; m=m'%输出m向量 运行结果: 图形如下: 三、求 的近似值,要求绝对误差小于 . 根据 解: 程序: (fuhua.m) functiony=fuhua(a,b)%输入为积分上下限,输出为积分值 n=1; t1=(b-a)/2*(exp(-a^2)/a+exp(-b^2)/b); while1 h=0; fori=0: n-1 h=exp(-((a+(2*i+1)*(b-a)/2/n)^ (2)))/(a+(2*i+1)*(b-a)/2/n)+h;%复化中矩形公式 end h2=(b-a)/n*h; y=1/(b-a)*(t1+h2); ifabs(y-t1)<0.5e-4;%判断绝对误差 break; end t1=y; n=2*n; end 运行结果: 四、构造迭代法求方程 的所有根,要求误差计算不超过 . 根据方程估算出两个解分别在0和1.5附近,然后进行计算。 解: 程序: (newton.m) x0=0.01;%估计第一个解在0附近 fori=1: 10000 x1=x0-(sin(x0)-(x0)^2/2)/(cos(x0)-x0); err=abs(x1-x0); x0=x1; if(err<10^(-10))%判断迭代误差 break; end end y1=x1%输出 x0=1.5;%估计第二个解在1.5附近 fori=1: NaN x1=x0-(sin(x0)-(x0)^2/2)/(cos(x0)-x0); err=abs(x1-x0); x0=x1; if(err<10^(-10))%判断迭代误差 break; end end y2=x1%输出 运行结果: 五、对于微分方程 用以下两种方法求解,步长 .画出数值解的图形并与精确解作比较. (1) 线性二步法 (2)经典Runge-Kutta法. 解: 将两种解法放在了一个程序里。 程序: (xianxing.m) formatlongg; a=0;b=1;h=0.05;%已知的初始值 fork=1: 1: 21 t(k)=a+(k-1)*h; true(k)=exp(-3*t(k))*11/9+2/3*t(k)-2/9;%根据解出的真值计算 end u (1)=1;%线性二步法 u (2)=u (1)+h*(2*t (1)-3*u (1)); fori=1: 1: 19 f(i)=2*t(i)-3*u(i); f(i+1)=2*t(i+1)-3*u(i+1); u(i+2)=u(i+1)+h*f(i+1); f(i+2)=2*t(i+2)-3*u(i+2); u(i+2)=u(i+1)+h/12*(5*f(i+2)+8*f(i+1)-f(i)); end u%输出线性二步法得出的值 urk (1)=1;%经典Runge-Kutta法 fori=1: 1: 20 k1=2*t(i)-3*u(i); k2=2*(t(i)+1/2*h)-3*(u(i)+1/2*h*k1); k3=2*(t(i)+1/2*h)-3*(u(i)+1/2*h*k2); k4=2*(t(i)+h)-3*(u(i)+h*k3); urk(i+1)=urk(i)+h/6*(k1+2*k2+2*k3+k4); end urk%输出经典R-K法得出的值 plot(t,true,'r')%用红色线画出真实值 holdon; plot(t,u,'k')%用黑色线画出线性二步法得出的值 holdon; plot(t,urk,'b-')%用蓝色点划线画出经典R-K法得出的值 fori=1: 1: 21 er1(i)=abs(u(i)-true(i));%计算线性二步法的误差 er2(i)=abs(urk(i)-true(i));%计算经典Runge-Kutta法的误差 end fori=1: 1: 21 er1(i)=u(i)-true(i);%计算线性二步法误差 er2(i)=urk(i)-true(i);%计算经典R-K法误差 end er1%输出线性二步法得出的误差值 er2%输出经典R-K法得出的误差值 运行结果: 两个方法的计算结果与真值比较表: N t(n) 精确解 U(tn) 线性二步法 经典Runge-Kutta法 数值解un 误差 数值解un 误差 0 0 1 1 0 1 0 1 0.05 0.86308752674 0.85 0.0130875267 0.8630882812 0.0000007545 2 0.1 0.74988893638 0.7395312500 .010******* 0.7517133203 0.0018243839 3 0.15 0.65710107420 0.6487905273 0.0083105468 0.6603687537 0.0032676795 4 0.2 0.58188088855 0.57525111846 0.0066297700 0.5863066369 0.0044257484 5 0.25 0.52178134223 0.5165256724 0.0052556696 0.5271309747 0.0053496325 6 0.3 0.47469625079 0.4705608179 0.0041354328 0.4807783093 0.0060820585 7 0.35 0.43881280446 0.4355877635 0.0032250408 0.4454712000 0.0066583955 8 0.4 0.41257070344 0.4100828975 0.0024878059 0.4196785835 0.0071078800 9 0.45 0.39462698523 0.3927338008 0.0018931843 0.4020816225 0.0074546373 10 0.5 0.38382575129 0.3824099738 0.0014157774 0.3915442885 0.0077185372 11 0.55 0.37917211053 0.3781376190 0.0010344914 0.3870880217 0.0079159112 12 0.6 0.37980975227 0.3790779208 0.0007318314 0.3878699041 0.0080601518 13 0.65 0.38500164305 0.3845083344 0.0004933085 0.3931638574 0.0081622144 14 0.7 0.39411341230 0.3938064697 0.0003069425 0.4023444477 0.0082310354 15 0.75 0.40659905224 0.4064362082 0.0001628439 0.4148729345 0.0082738822 16 0.8 0.42198860957 0.4219357447 0.0000538648 0.4302852541 0.0082966445 17 0.85 0.43987759177 0.4399072864 0.0000296946 0.4481816684 0.0082874677 18 0.9 0.45991784890 0.4600081800 0.0000903310 0.4682178482 0.0083040766 19 0.95 0.48180972551 0.4819432685 0.0001335430 0.4900971932 0.0082999993 20 1.0 0.50529530578 0.5054583077 0.0001630019 0.5135642157 0.0082689099
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 矩阵 报告