矩阵与数值分析例题matlab仿真
- 格式:pdf
- 大小:414.19 KB
- 文档页数:9
数值分析中求解线性方程组的MATLAB程序(6种)1.回溯法(系数矩阵为上三角)function X=uptrbk(A,B)%求解方程组,首先化为上三角,再调用函数求解[N,N]=size(A);X=zeros(N,1);C=zeros(1,N+1);Aug=[A B];for p=1:N-1[Y,j]=max(abs(Aug(p:N,p)));C=Aug(p,:);Aug(p,:)=Aug(j+p-1,:);Aug(j+p-1,:)=C;if Aug(p,p)==0'A was singular.No unique solution.'break;endfor k=p+1:Nm=Aug(k,p)/Aug(p,p);Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);endendD=Aug;X=backsub(Aug(1:N,1:N),Aug(1:N,N+1));2.系数矩阵为下三角function x=matrix_down(A,b)%求解系数矩阵是下三角的方程组n=length(b);x=zeros(n,1);x(1)=b(1)/A(1,1);for k=2:1:nx(k)=(b(k)-A(k,1:k-1)*x(1:k-1))/A(k,k);end3.普通系数矩阵(先化为上三角,在用回溯法)function X=uptrbk(A,B)%求解方程组,首先化为上三角,再调用函数求解[N,N]=size(A);X=zeros(N,1);C=zeros(1,N+1);Aug=[A B];for p=1:N-1[Y,j]=max(abs(Aug(p:N,p)));C=Aug(p,:);Aug(p,:)=Aug(j+p-1,:);Aug(j+p-1,:)=C;if Aug(p,p)==0'A was singular.No unique solution.'break;endfor k=p+1:Nm=Aug(k,p)/Aug(p,p);Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);endendD=Aug;X=backsub(Aug(1:N,1:N),Aug(1:N,N+1));4.三角分解法function [X,L,U]=LU_matrix(A,B)%A是非奇异矩阵%AX=B化为LUX=B,L为下三角,U为上三角%程序中并没有真正解出L和U,全部存放在A中[N,N]=size(A);X=zeros(N,1);Y=zeros(N,1);C=zeros(1,N);R=1:N;for p=1:N-1[max1,j]=max(abs(A(p:N,p)));C=A(p,:);A(p,:)=A(j+p-1,:);A(j+p-1,:)=C;d=R(p);R(p)=R(j+p-1);R(j+p-1)=d;if A(p,p)==0'A is singular.No unique solution'break;endfor k=p+1:Nmult=A(k,p)/A(p,p);A(k,p)=mult;A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N);endendY(1)=B(R(1));for k=2:NY(k)=B(R(k))-A(k,1:k-1)*Y(1:k-1);endX(N)=Y(N)/A(N,N);for k=N-1:-1:1X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k);endL=tril(A,-1)+eye(N)U=triu(A)5.雅克比迭代法function X=jacobi(A,B,P,delta,max1);%雅克比迭代求解方程组N=length(B);for k=1:max1for j=1:NX(j)=(B(j)-A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N]))/A(j,j);enderr=abs(norm(X'-P));relerr=err/(norm(X)+eps);P=X';if (err<delta)|(relerr<delta)breakendendX=X';k6.盖斯迭代法function X=gseid(A,B,P,delta,max1);%盖斯算法,求解赋初值的微分方程N=length(B);for k=1:max1for j=1:Nif j==1X(1)=(B(1)-A(1,2:N)*P(2:N))/A(1,1);elseif j==NX(N)=(B(N)-A(N,1:N-1)*(X(1:N-1))')/A(N,N);elseX(j)=(B(j)-A(j,1:j-1)*X(1:j-1)-A(j,j+1:N)*P(j+1:N))/A(j,j);endenderr=abs(norm(X'-P));relerr=err/(norm(X)+eps);P=X';if (err<delta)|(relerr<delta)break;endendX=X';k。
(Hilbert 矩阵)病态线性方程组的求解理论分析表明,数值求解病态线性方程组很困难。
考虑求解如下的线性方程组的求解Hx = b ,期中H 是Hilbert 矩阵,()ij n n Hh ,11ij h i j ,i ,j = 1,2,…,n 1.估计矩阵的2条件数和阶数的关系2.对不同的n ,取(1,1,,1)nx K ?,分别用Gauss 消去,Jacobi 迭代,Gauss-seidel 迭代,SOR 迭代和共轭梯度法求解,比较结果。
3.结合计算结果,试讨论病态线性方程组的求解。
第1小题:condition.m %第1小题程序t1=20;%阶数n=20x1=1:t1;y1=1:t1;for i=1:t1H=hilb(i);y1(i)=log(cond(H));endplot(x1,y1);xlabel('阶数n');ylabel('2-条件数的对数(log(cond(H))');title('2-条件数的对数(log(cond(H))与阶数n 的关系图');t2=200;%阶数n=200x2=1:t2;y2=1:t2;for i=1:t2H=hilb(i);y2(i)=log(cond(H));endplot(x2,y2);xlabel('阶数n');ylabel('2-条件数的对数(log(cond(H))');title('2-条件数的对数(log(cond(H))与阶数n 的关系图');画出Hilbert 矩阵2-条件数的对数和阶数的关系n=200时n=20时从图中可以看出,1)在n小于等于13之前,图像近似直线log(cond(H))~1.519n-1.8332)在n大于13之后,图像趋于平缓,并在一定范围内上下波动,同时随着n的增加稍有上升的趋势第2小题:solve.m%m第2小题主程序N=4000;xGauss=zeros(N,1);xJacobi=zeros(N,1);xnJ=zeros(N,1);xGS=zeros(N,1);xnGS=zeros(N,1);xSOR=zeros(N,1);xnSOR=zeros(N,1);xCG=zeros(N,1);xnCG=zeros(N,1);for n=1:N;x=ones(n,1);t=1.1;%初始值偏差x0=t*x;%迭代初始值e=1.0e-8;%给定的误差A=hilb(n);b=A*x;max=100000000000;%可能最大的迭代次数w=0.5;%SOR迭代的松弛因子G=Gauss(A,b);[J,nJ]=Jacobi(A,b,x0,e,max);[GS,nGS]=G_S(A,b,x0,e,max);[S_R,nS_R]=SOR(A,b,x0,e,max,w);[C_G,nC_G]=CG(A,b,x0,e,max);normG=norm(G'-x);xGauss(n)=normG;normJ=norm(J-x);nJ;xJacobi(n)=normJ;xnJ(n)=nJ;normGS=norm(GS-x);nGS;xGS(n)=normGS;xnGS(n)=nGS;normS_R=norm(S_R-x);nS_R;xSOR(n)=normS_R;xnSOR(n)=nS_R;normC_G=norm(C_G-x);nC_G;xCG(n)=normC_G;xnCG(n)=nC_G;endGauss.m%Gauss消去法function x=Gauss(A,b)n=length(b);l=zeros(n,n);x=zeros(1,n);%消去过程for i=1:n-1for j=i+1:nl(j,i)=A(j,i)/A(i,i);for k=i:nA(j,k)=A(j,k)-l(j,i)*A(i,k);endb(j)=b(j)-l(j,i)*b(i);endend%回代过程x(n)=b(n)/A(n,n);for i=n-1:-1:1c=A(i,:).*x;x(i)=(b(i)-sum(c(i+1:n)))/A(i,i);endJacobi.m%Jacobi迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m 可能最大的迭代次数function [x,n]=Jacobi(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=D\(L+U);f=D\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Jacobi迭代次数过多,迭代可能不收敛');break;endendG_S.m%Gauss-Seidel迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=G_S(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-L)\U;f=(D-L)\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Gauss-Seidel迭代次数过多,迭代可能不收敛');break;endendSOR.m%SOR超松弛迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数,w松弛因子function [x,n]=SOR(A,b,x0,e,m,w)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-w*L)\((1-w)*D+w*U);f=(D-w*L)\b*w;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('SOR超松弛迭代次数过多,迭代可能不收敛');break;endendCG.m%CG共轭梯度法,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=CG(A,b,x0,e,m)r=b-A*x0;p=r;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=1;while norm(r1)>ebelta=(r1'*r1)/(r'*r);p=r1+belta*p;r=r1;x0=x;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=n+1;if n>mdisp('CG共轭梯度法迭代次数过多,迭代可能不收敛');break;endend。
1.1 矩阵的表示1.1.1 数值矩阵的生成1.实数值矩阵输入MATLAB的强大功能之一体现在能直接处理向量或矩阵。
当然首要任务是输入待处理的向量或矩阵。
不管是任何矩阵(向量),我们可以直接按行方式输入每个元素:同一行中的元素用逗号(,)或者用空格符来分隔,且空格个数不限;不同的行用分号(;)分隔。
所有元素处于一方括号([])内;当矩阵是多维(三维以上),且方括号内的元素是维数较低的矩阵时,会有多重的方括号。
如:>> Time = [11 12 1 2 3 4 5 6 7 8 9 10]Time =11 12 1 2 3 4 5 6 7 8 9 10>> X_Data = [2.32 3.43;4.37 5.98]X_Data =2.433.434.375.98>> vect_a = [1 2 3 4 5]vect_a =1 2 3 4 5>> Matrix_B = [1 2 3;>> 2 3 4;3 4 5]Matrix_B = 1 2 32 3 43 4 5>> Null_M = [ ] %生成一个空矩阵2.复数矩阵输入复数矩阵有两种生成方式:第一种方式例1-1>> a="2".7;b=13/25;>> C=[1,2*a+i*b,b*sqrt(a); sin(pi/4),a+5*b,3.5+1]C=1.0000 5.4000 + 0.5200i 0.85440.7071 5.3000 4.5000第2种方式例1-2>> R=[1 2 3;4 5 6], M=[11 12 13;14 15 16]R =1 2 34 5 6M =11 12 1314 15 16>> CN="R"+i*MCN =1.0000 +11.0000i2.0000 +12.0000i3.0000 +13.0000i4.0000 +14.0000i5.0000 +15.0000i6.0000 +16.0000i1.1.2 符号矩阵的生成在MATLAB中输入符号向量或者矩阵的方法和输入数值类型的向量或者矩阵在形式上很相像,只不过要用到符号矩阵定义函数sym,或者是用到符号定义函数syms,先定义一些必要的符号变量,再像定义普通矩阵一样输入符号矩阵。
《矩阵与数值分析》上机大作业1.给定n 阶方程组Ax b =,其中6186186186A ⎛⎫ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪⎝⎭ ,7151514b ⎛⎫ ⎪ ⎪ ⎪= ⎪⎪⎪⎝⎭则方程组有解(1,1,,1)Tx = 。
对10n =和84n =,分别用Gauss 消去法和列主元消去法解方程组,并比较计算结果。
%产生三对角矩阵 n=84; %或n=10;A=zeros(n); b=zeros(1,n); for i=1:n-1A(i,i)=6;A(i,i+1)=1;A(i+1,i)=8; endA(n,n)=6;for i=2:n-1 b(1)=6; b(i)=15; b(n)=14; end Ab=[A b'];%Gauss 消元法for j=1:n-1 %按列循环 for k=j+1:n %消元Ab(k,:)=Ab(k,:)-Ab(j,:)*(Ab(k,j)/Ab(j,j)); end endx(n)=Ab(n,n+1)/Ab(n,n); for i=n-1:-1:1 %回代法求x for j=n:-1:i+1Ab(i,n+1)=Ab(i,n+1)-Ab(i,j)*x(j); endx(i)=Ab(i,n+1)/Ab(i,i); end(1)当n=10时,Gauss 消去法 Gauss 列主元消去法 x=1.000000000000000 x=1.000000000000000 1.000000000000000 1.0000000000000001.000000000000000 1.000000000000000 1.000000000000001 1.000000000000000 0.999999999999998 1.000000000000000 1.000000000000004 1.000000000000000 0.999999999999993 1.000000000000000 1.000000000000012 1.000000000000000 0.999999999999979 1.000000000000000 1.000000000000028 1.000000000000000(2) 当n=84时,Gauss 消去法的解是错解Columns 34 through 392147483649.00000 -4294967295.00000 8589934592.99999 -17179869182.9999 34359738368.9998Gauss 列主元消去法x 与x=(1,1…1)T 偏差不大 Columns 34 through 391.000000172108412 0.999999661246936 1.000000655651093 0.999998776117961 1.000002098083496综上,高斯列主元消去法可以避免小数作除数带来的误差,获得满意的数值解。
2.12 MATLAB 语言的数值运算-实训2.12.1实训目的1.学会矩阵的建立方法及其矩阵的转置、相乘、求逆等运算;2.识别了解特殊矩阵;3.学会求解方程与方程组;4.学会通过编程解决一些实际问题; 2.12.2实训内容 1.矩阵建立及其运算]4321012345[1-----=A[]0.19.08.07.06.05.04.03.02.01.02=A ]0.19.08.07.06.05.04.03.02.01.0[3=A][40.19.08.07.06.05.04.03.02.01.0e e e e e e e e e e A = ]2222222222[50.19.08.07.06.05.04.03.02.01.0=A求(1)211A A D += (2)232A A D -= (3)4/.13A D = (4)5*.44A A D = (5)5*35A D = (6))2(.^16A D =>> A1=-5:4;>> A2=0.1:0.1:1.0; >> A3=sqrt(A2); >> A4=exp(A2); >> A5=2.^(A2); >> format bank >> D1=A1+A2 D1 =-4.90 -3.80 -2.70 -1.60 -0.50 0.60 1.70 2.80 3.90 5.00 >> D2=A3-A2 D2 =0.22 0.25 0.25 0.23 0.21 0.17 0.14 0.09 0.05 0 >> D3=1./A4 D3 =0.90 0.82 0.74 0.67 0.61 0.55 0.50 0.45 0.41 0.37 >> D4=A4.*A5D4 =1.18 1.40 1.66 1.972.33 2.763.27 3.874.595.44 >> D5=3*A5 D5 =3.22 3.45 3.69 3.964.24 4.55 4.875.22 5.606.00 >> D6=A1.^2 D6 =25.00 16.00 9.00 4.00 1.00 0 1.00 4.00 9.00 16.002.建立矩阵⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-----=1418712106114231359152345686420B(1)矩阵B 的逆矩阵)(inv B (2)矩阵B 对应的行列式)det(B>> B=[0:2:8;-6:-2;15,9,5,13,3;2,4,11,6,10;12,7,8,1,14]B =0 2.00 4.00 6.00 8.00 -6.00 -5.00 -4.00 -3.00 -2.00 15.00 9.00 5.00 13.00 3.00 2.00 4.00 11.00 6.00 10.00 12.00 7.00 8.00 1.00 14.00 >> inv(B) ans =-0.18 0.29 0.11 0.07 0.07 0.31 -0.72 -0.20 -0.23 -0.07 -0.21 0.11 0.03 0.21 -0.02 0.05 0.11 0.08 0.02 -0.04 0.12 0.04 -0.02 -0.06 0.06 >> det(B) ans =-17568.00(3)利用矩阵元素的提取方法建立以下矩阵 ①矩阵b01:矩阵B 的3~4行元素; ②矩阵b02:矩阵B 的2~5列元素;③矩阵b03:由矩阵B 的1~3行2~4列交叉点所对应的元素组成; >> B=[0:2:8;-6:-2;15,9,5,13,3;2,4,11,6,10;12,7,8,1,14]; >> b01=B(3:4,:) b01 =15.00 9.00 5.00 13.00 3.002.00 4.00 11.00 6.00 10.00 >> b02=B(:,2:5) b02 =2.00 4.00 6.00 8.00 -5.00 -4.00 -3.00 -2.00 9.00 5.00 13.00 3.004.00 11.00 6.00 10.00 7.00 8.00 1.00 14.00 >> b03=B(1:3,2:4) b03 =2.00 4.00 6.00 -5.00 -4.00 -3.00 9.00 5.00 13.00 3.矩阵的转置与翻转已知矩阵⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=151413121110987654321m ,求取以下矩阵观察并记录。
matlab仿真实例100题Matlab是一种强大的数学软件,广泛应用于科学计算、数据分析和工程仿真等领域。
在学习和使用Matlab的过程中,通过实例的方式进行仿真练习是一种非常有效的学习方法。
下面将给出100个Matlab仿真实例题目,帮助读者更好地掌握Matlab的使用。
1. 编写一个程序,计算并输出1到100之间所有奇数的和。
2. 编写一个程序,计算并输出1到100之间所有偶数的乘积。
3. 编写一个程序,计算并输出1到100之间所有素数的个数。
4. 编写一个程序,计算并输出1到100之间所有整数的平方和。
5. 编写一个程序,计算并输出1到100之间所有整数的立方和。
6. 编写一个程序,计算并输出1到100之间所有整数的阶乘和。
7. 编写一个程序,计算并输出1到100之间所有整数的倒数和。
8. 编写一个程序,计算并输出1到100之间所有整数的平均值。
9. 编写一个程序,计算并输出1到100之间所有整数的中位数。
10. 编写一个程序,计算并输出1到100之间所有整数的标准差。
11. 编写一个程序,计算并输出1到100之间所有整数的方差。
12. 编写一个程序,计算并输出1到100之间所有整数的最大值。
13. 编写一个程序,计算并输出1到100之间所有整数的最小值。
15. 编写一个程序,计算并输出1到100之间所有整数的平方根和。
16. 编写一个程序,计算并输出1到100之间所有整数的立方根和。
17. 编写一个程序,计算并输出1到100之间所有整数的对数和。
18. 编写一个程序,计算并输出1到100之间所有整数的指数和。
19. 编写一个程序,计算并输出1到100之间所有整数的正弦和。
20. 编写一个程序,计算并输出1到100之间所有整数的余弦和。
21. 编写一个程序,计算并输出1到100之间所有整数的正切和。
22. 编写一个程序,计算并输出1到100之间所有整数的双曲正弦和。
23. 编写一个程序,计算并输出1到100之间所有整数的双曲余弦和。
一、最小二乘法,用MATLAB实现1. 数值实例下面给定的是乌鲁木齐最近1个月早晨7:00左右(新疆时间)的天气预报所得到的温度,按照数据找出任意次曲线拟合方程和它的图像。
下面用MATLAB编程对上述数据进行最小二乘拟合。
下面用MATLAB编程对上述数据进行最小二乘拟合2、程序代码x=[1:1:30];y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7,6,5,3,1];a1=polyfit(x,y,3) %三次多项式拟合%a2= polyfit(x,y,9) %九次多项式拟合%a3= polyfit(x,y,15) %十五次多项式拟合%b1=polyval(a1,x)b2=polyval(a2,x)b3=polyval(a3,x)r1= sum((y-b1).^2) %三次多项式误差平方和%r2= sum((y-b2).^2) %九次次多项式误差平方和%r3= sum((y-b3).^2) %十五次多项式误差平方和%plot(x,y,'*') %用*画出x,y图像%hold onplot(x,b1, 'r') %用红色线画出x,b1图像%hold onplot(x,b2, 'g') %用绿色线画出x,b2图像%hold onplot(x,b3, 'b:o') %用蓝色o线画出x,b3图像%3、数值结果不同次数多项式拟合误差平方和为:r1=67.6659r2=20.1060r3=3.7952r1、r2、r3分别表示三次、九次、十五次多项式误差平方和。
4、拟合曲线如下图二、 线性方程组的求解( 高斯-塞德尔迭代算法 )1、实例: 求解线性方程组(见书P233页)⎪⎪⎩⎪⎪⎨⎧=++=-+=+-3612363311420238321321321x x x x x x x x x 记A x=b, 其中⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=363320,,12361114238321b x A x x x任取初始值()()Tx0000=,进行迭代。
利用Matlab解决常见数学问题的案例分析概述:Matlab是一款流行的科学软件,广泛应用于数学建模、数据分析、图像处理等领域。
本文将通过几个实际案例,介绍如何利用Matlab解决常见的数学问题,并分析其解决方法和效果。
案例一:线性方程组的求解线性方程组是数学中常见的问题之一。
假设有如下线性方程组:3x + 2y = 14x - 3y = 5可以使用Matlab中的线性方程组求解函数`linsolve`来求解。
首先,定义系数矩阵A和常数矩阵b,并调用`linsolve`函数求解方程组:```matlabA = [3 2; 4 -3];b = [1; 5];x = linsolve(A, b);```运行上述代码后,可以得到方程组的解x为:x = 3y = -2案例二:函数曲线绘制Matlab具有强大的绘图功能,可以绘制各种函数曲线。
例如,我们可以绘制正弦函数sin(x)在区间[-2π,2π]上的曲线。
首先,定义x的取值范围,并计算对应的y 值:```matlabx = -2*pi:0.1:2*pi;y = sin(x);```接下来,使用`plot`函数将曲线绘制出来:```matlabplot(x, y);```运行代码后,可以得到正弦函数的曲线图。
案例三:最小二乘拟合最小二乘拟合是一种常见的曲线拟合方法,用于将一组数据拟合成一条曲线。
假设有一组离散的数据点,我们希望找到一个曲线来拟合这些数据。
在Matlab中,可以使用`polyfit`函数进行最小二乘拟合。
例如,假设有一组数据:x = [1 2 3 4 5];y = [0.5 2.5 2 4 3.5];可以使用`polyfit`函数进行线性拟合:```matlabp = polyfit(x, y, 1);```其中,第一个参数x是自变量的取值,第二个参数y是因变量的取值,第三个参数1表示进行一次多项式拟合。
拟合的结果保存在向量p中,p(1)为拟合曲线的斜率,p(2)为截距。
实验二 MATLAB 矩阵分析与处理班级 姓名 学号1. 设有分块矩阵33322322E R A O S ⨯⨯⨯⨯⎡⎤=⎢⎥⎣⎦,其中E 、R 、O 、S 分别为单位矩阵、随机矩阵、零矩阵和对角矩阵,试通过数值计算验证22+E R R S A O S ⎡⎤=⎢⎥⎣⎦。
>> E=eye(3); >> R=rand(3,2); >> O=zeros(2,3); >> S=diag([1 2]); >> A=[E R;O S]; >> A^2>> B=[E R+R*S;O S^2]2. 产生5阶希尔伯特矩阵H 和5阶帕斯卡矩阵P ,且求其行列式的值Hh 和Hp 以及它们的条件数Th 和Tp ,判断哪个矩阵性能更好。
为什么? >> H=hilb(5); >> P=pascal(5); >> Hp=det(P); >> Hh=det(H);>> Th=cond(H,2) >> Tp=cond(P ,2)3. 建立一个55⨯矩阵,求它的行列式值、迹、秩和范数。
>> A=rand(5); >> det(A); >> trace(A); >> rank(A); >> norm(A,1); >> norm(A,2); >> norm(A,inf);4. 已知2961820512885A -⎡⎤⎢⎥=⎢⎥⎢⎥-⎣⎦求A 的特征值及特征向量。
>> A=[-29 6 18;20 5 12;-8 8 5]; >> [V ,D]=eig(A) V =0.7130 0.2803 0.2733 -0.6084 -0.7867 0.8725 0.3487 0.5501 0.4050 D =-25.3169 0 0 0 -10.5182 0 0 0 16.83515. 下面是一个线性方程组1231/21/31/40.951/31/41/50.671/41/51/60.52x x x ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦(1)求方程的解。
1、 在MATLAB 中用Jacobi 迭代法讨论线性方程组,1231231234748212515x x x x x x x x x -+=⎧⎪-+=-⎨⎪-++=⎩(1)给出Jacobi 迭代法的迭代方程,并判定Jacobi 迭代法求解此方程组是否收敛。
(2)若收敛,编程求解该线性方程组.解(1):A=[4 -1 1;4 —8 1;-2 1 5] %线性方程组系数矩阵A =4 -1 1 4 -8 1 —2 1 5>> D=diag(diag(A))D =4 0 0 0 —8 0 0 0 5〉〉 L=—tril (A,-1) % A 的下三角矩阵L =0 0 0 —4 0 0 2 —1 0〉〉U=-triu(A,1)% A的上三角矩阵U =0 1 —10 0 —10 0 0B=inv(D)*(L+U)% B为雅可比迭代矩阵B =0 0.2500 —0。
25000.5000 0 0.12500。
4000 —0.2000 0〉〉r=eigs(B,1)%B的谱半径r =0。
3347 〈1Jacobi迭代法收敛。
(2)在matlab上编写程序如下:A=[4 —1 1;4 -8 1;—2 1 5];〉〉b=[7 —21 15]';>〉x0=[0 0 0]’;〉〉[x,k]=jacobi(A,b,x0,1e—7)x =2。
00004.00003。
0000k =17附jacobi迭代法的matlab程序如下:function [x,k]=jacobi(A,b,x0,eps)% 采用Jacobi迭代法求Ax=b的解%A为系数矩阵%b为常数向量%x0为迭代初始向量%eps为解的精度控制max1= 300; %默认最多迭代300,超过300次给出警告D=diag(diag(A));%求A的对角矩阵L=-tril(A,—1); %求A的下三角阵U=—triu(A,1); %求A的上三角阵B=D\(L+U);f=D\b;x=B*x0+f;k=1;%迭代次数while norm(x-x0)>=epsx0=x;x=B*x0+f;k=k+1;if(k〉=max1)disp(’迭代超过300次,方程组可能不收敛’);return;endend2、设有某实验数据如下:(1)在MATLAB中作图观察离散点的结构,用多项式拟合的方法拟合一个合适的多项式函数;(2)在MATLAB中作出离散点和拟合曲线图。
2013级工科硕士研究生《矩阵与数值分析》课程数值实验题目一、设622101NNjSj==-∑,分别编制从小到大和从大到小的顺序程序分别计算100001000000,S S并指出两种方法计算结果的有效位数。
Matlab程序如下:function [si,sd]=S(N)format long;si=0;sd=0;for j=N:-1:2si=1.0e6/(j^2-1)+si;endfor j=2:Nsd=1.0e6/(j^2-1)+sd;endend在matlab命令窗口中输入:[si,sd]=S(10000)运行结果:si =7.499000049995000e+005sd =7.499000049994994e+005在matlab命令窗口中输入:[si,sd]=S(1000000)运行结果:si =7.499990000005000e+005sd =7.499990000005200e+0051结果分析:si为从大到小的顺序求和的值,sd为从小到大的顺序求和的值。
当N分别为10000和1000000时,si分别为7.499000049995000e+005和7.499990000005000e+005,可以看出这两个数的有效值均为13位;而sd分别为7.499000049994994e+005和7.499990000005200e+005,这两个数的有效值均为16位。
这就出现了我们在矩阵理论课上所学的“大数吃小数”的问题。
为了使结果更为精确我们必须避免在四则运算中出现“大数吃小数”的情况,应该按从小到大的顺序进行求和。
二、解线性方程组1.分别利用Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组Ax b=,其中常向量为()21n-维随机生成的列向量,系数矩阵A具有如下形式1111111122n n n n n n n n T I I I A I I T I --------+-⎛⎫ ⎪- ⎪= ⎪- ⎪-+⎝⎭O O O O , 其中1211112n T --⎛⎫ ⎪- ⎪= ⎪- ⎪-⎝⎭O O O O 为1n -阶矩阵,1n I -为1n -阶单位矩阵,迭代法计算停止的条件为:101210k k x x -+-<,给出10,100,1000n =时的不同迭代步数. 求解系数矩阵A 和向量b 的Matlab 程序如下:function [A,b,x0]=jz(n)for i=1:n-1 %此处n 赋值即可,如n=100for j=1:n-1if(i==j)T(i,j)=2;endif(i==(j+1))T(i,j)=-1;endif(i==(j-1))T(i,j)=-1;endendendI=eye(n-1);k=1;for mm=1:(n-1)A(k:(k+n-2),k:(k+n-2))=T+2*I;k=k+n-1;endk=1;for xx=1:(n-2)A(k:(k+n-2),(k+n-1):(k+2*n-3))=-I;A((k+n-1):(k+2*n-3),k:(k+n-2))=-I;k=k+n-1;endb=randn((n-1)^2,1);x0=zeros((n-1)^2,1);Jacobi 迭代法Matlab 程序如下:function n=jacobi(A,b,x0)eps= 1.0e-10;M = 10000;D=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵B=D\(L+U);f=D\b;x=B*x0+f;n=1; %迭代次数while norm(x-x0)>=epsx0=x;x=B*x0+f;n=n+1;if(n>=M)disp('Warning: 迭代次数太多,可能不收敛!');return;endendGauss-Seidel迭代法Matlab程序如下:function n=gauseidel(A,b,x0)eps= 1.0e-10;M = 10000;D=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵G=(D-L)\U;f=(D-L)\b;x=G*x0+f;n=1; %迭代次数while norm(x-x0)>=epsx0=x;x=G*x0+f;n=n+1;if(n>=M)disp('Warning: 迭代次数太多,可能不收敛!');return;endend在matlab命令窗口中输入:[A,b,x0]=jz(10);J10=jacobi(A,b,x0)G10=gauseidel(A,b,x0)[A,b,x0]=jz(20);J20=jacobi(A,b,x0)G20=gauseidel(A,b,x0)[A,b,x0]=jz(30);J30=jacobi(A,b,x0)G30=gauseidel(A,b,x0)运行结果:J10 =436;G10 =214J20 =1810;G20 =913J30 =3990;G30 =2001结果分析:N=10且M=1000时,Jacobi 迭代法和Gauss —seidel 迭代法的迭代次数分别为436和214;N=20且M=1000时,Jacobi 迭代法和Gauss —seidel 迭代法的迭代次数分别为1810和913;N=30且M=1000时 ,Jacobi 和Gauss-seidel 算法的迭代次数分别为3990和2001次。
数值分析实验报告学院:电气工程与自动化学院专业:控制理论与控制工程姓名:李亚学号:61201401622014 年 12 月24日实验一 函数插值方法一、目的和意义1、 学会常用的插值方法,求函数的近似表达式,以解决其它实际问题;2、 明确插值多项式和分段插值多项式各自的优缺点;3、 熟悉插值方法的程序编制;4、 如果绘出插值函数的曲线,观察其光滑性。
二、实验原理1、 Lagrange 插值公式00,()n ni n k k i i k k i x x L x y x x ==≠⎛⎫-= ⎪-⎝⎭∑∏编写出插值多项式程序;2、 给出插值多项式或分段三次插值多项式的表达式;三、实验要求对于给定的一元函数)(x f y =的n+1个节点值(),0,1,,j j y f x j n ==。
试用Lagrange 公式求其插值多项式或分段二次Lagrange 插值多项式。
数据如下:(1求五次Lagrange 多项式5L ()x ,计算(0.596)f ,(0.99)f 的值。
(提示:结果为(0.596)0.625732f ≈, (0.99) 1.05423f ≈)试构造Lagrange 多项式6,和分段三次插值多项式,计算的(1.8)f ,(6.15)f 值。
(提示:结果为(1.8)0.164762f ≈, (6.15)0.001266f ≈)四、实验过程1.进入matlab开发环境;2.根据实验内容和要求编写程序,程序如下所示,程序通过运用function 函数编写,生成.m文件。
调用时只需要在命令窗口调用y=Lagrange(A,input)就可以实现任意次数拉格朗日插值法求解。
function y=Lagrange(A,input)[a,b]=size(A);x=input;y=0;for j=1:aMj=1;Nj=1;for k=1:aif(k==j)continue;endMj=Mj*(x-A(k,1));Nj=Nj*(A(j,1)-A(k,1));endy=y+A(j,2)*Mj/Nj;end3.调试程序并运行程序;调用拉格朗日脚本文件对以上两个表格数据求解,表格一对应MATLAB向量A;表格二对应向量I。
标题:深度解析:如何使用MATLAB解矩阵方程组的例题在数学和工程领域中,矩阵方程组是一种常见的问题,解决这类问题的方法有很多种,其中MATLAB提供了方便而强大的工具。
本文将深入探讨如何使用MATLAB来解决矩阵方程组的例题,以便读者能更深入地理解和掌握这一重要的数学概念。
矩阵方程组是指包含矩阵和向量的一组线性方程,一般形式为Ax=b,其中A是系数矩阵,x是未知向量,b是已知向量。
解矩阵方程组的核心问题是找到满足方程的x,而MATLAB提供了多种方法来解这类问题。
接下来,我们将通过一个具体的例题来详细介绍MATLAB的解决方法。
假设有如下的矩阵方程组:```2x + 3y - z = 14x - 2y + 3z = -23x + y + 2z = 3```我们可以将系数矩阵A和常数向量b定义为MATLAB中的矩阵和向量:```A = [2, 3, -1; 4, -2, 3; 3, 1, 2];b = [1; -2; 3];```接下来,我们可以使用MATLAB中的反斜杠运算符`\`来求解矩阵方程组,即x=A\b。
在MATLAB中,这个简单的运算就能够求解出矩阵方程组的解向量x。
通过这种方式,我们可以非常方便地求解包含大量线性方程的矩阵方程组,而无需手动进行繁琐的计算。
在MATLAB中,还有一些其他函数和工具可以用来解决矩阵方程组,如inv函数和linsolve函数等。
这些工具在不同情况下可能会有不同的性能表现和数值稳定性,因此在实际应用中需要根据具体情况进行选择和使用。
总结来说,MATLAB提供了多种简单而强大的工具来解决矩阵方程组,无论是小规模的方程组还是大规模的方程组,都能够得到高效和准确的解。
通过本文的详细介绍和例题分析,相信读者已经对MATLAB解决矩阵方程组有了更深入的理解和掌握。
在解决矩阵方程组的过程中,我个人认为MATLAB的丰富函数库和简洁的语法是其最大的优势之一。
无论是初学者还是专业人士,都可以借助MATLAB轻松地解决复杂的矩阵运算问题。
1.使用444(x)对数据进行插值,并写出误差分析理论。
建立脚本x1=[0.2 0.4 0.6 0.8 1.0];y1=[0.98 0.92 0.81 0.64 0.38];n=length(y1);c=y1(:);for j=2:n %求差商for i=n:-1:jc(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1));endendsyms x df d;df(1)=1;d(1)=y1(1);for i=2:n %求牛顿差值多项式df(i)=df(i-1)*(x-x1(i-1));d(i)=c(i)*df(i);enddisp('4次牛顿插值多项式');P4=vpa(collect((sum(d))),5) %P4即为4次牛顿插值多项式,保留小数点后5位数figureezplot(P4,[0.2,1.08]);输出结果为P4 =- 0.52083*x^4 + 0.83333*x^3 - 1.1042*x^2 + 0.19167*x + 0.98插值余项:R4(x)=f(5)( ξ)/ (5!)* (x - 0.6)*(x - 0.4)*(x - 0.8)*(x - 1)*(x-0.2)新建一个M-filesyms x l;x1=[0.2 0.4 0.6 0.8 1.0];y1=[0.98 0.92 0.81 0.64 0.38];n=length(x1);Ls=sym(0);for i=1:nl=sym(y1(i));for k=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfor k=i+1:nl=l*(x-x1(k))/(x1(i)-x1(k));endLs=Ls+l;endLs=simplify(Ls) %为所求插值多项式Ls(x).figureezplot(Ls,[0.2,1.08]);输出结果为Ls =- (25*x^4)/48 + (5*x^3)/6 - (53*x^2)/48 + (23*x)/120 + 49/50插值余项:R4(x)=f(5)( ξ)/ (5!)* (x - 0.6)*(x - 0.4)*(x - 0.8)*(x - 1)*(x-0.2) 2.试求3次、合曲线,用图示数据曲线及相应的三种拟合曲线。