当前位置：文档之家 › 矩阵与数值分析例题matlab仿真

矩阵与数值分析例题matlab仿真

  • 格式：pdf
  • 大小：414.19 KB
  • 文档页数：9

下载文档原格式

  / 9

合集下载

数值分析中求解线性方程组的MATLAB程序(6种)

数值分析中求解线性方程组的MATLAB程序(6种)

数值分析中求解线性方程组的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。

matlab在数值分析中的应用5

matlab在数值分析中的应用5

• ym(P) 或 f=poly2sym(P,x)
格式：
P=sym2poly(f)
• 例：
>> P=[1 2 3 4 5 6]; % 先由系数按降幂顺序排列表示多 项式 >> f=poly2sym(P,'v') % 以 v 为算子表示多项式 f= v^5+2*v^4+3*v^3+4*v^2+5*v+6
>> P=sym2poly(f) P= 1 2 3 4 5
6
• 矩阵的逆矩阵
格式：
C=inv(A)
例：
>> format long; H=hilb(4); H1=inv(H)
H1 = 1.0e+003 * 0.01600000000000 -0.11999999999999 0.23999999999998 -0.13999999999999
>> H=sym(hilb(30)); norm(double(H*inv(H)-eye(size(H))))
ans = 0
• 例：奇异阵求逆
>> A=[16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1]; >> format long; B = inv(A)
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.306145e-017. B= 1.0e+014 * 0.93824992236885 2.81474976710656 -2.81474976710656 2.81474976710656 8.44424930131968 -8.44424930131968 -2.81474976710656 -8.44424930131968 8.44424930131968 -0.93824992236885 -2.81474976710656 2.81474976710656

数值分析(hilbert矩阵)病态线性方程组的求解matlab程序

数值分析(hilbert矩阵)病态线性方程组的求解matlab程序

（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。

matlab 矩阵相关操作及实例

matlab 矩阵相关操作及实例

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，先定义一些必要的符号变量，再像定义普通矩阵一样输入符号矩阵。

《矩阵与数值分析》上机大作业matlab

《矩阵与数值分析》上机大作业matlab

《矩阵与数值分析》上机大作业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综上，高斯列主元消去法可以避免小数作除数带来的误差，获得满意的数值解。

MATLAB语言及控制系统仿真_参考答案解析_第2章

MATLAB语言及控制系统仿真_参考答案解析_第2章

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上机作业

x=zeros(n,1); %为列向量 x 预分配存储空间 y=1:n; %定义行向量 y y=y'; %把行向量 y 改成列向量 for i=1:n
x(i)=1/i; %按要求给向量 x 赋值，其值递减 end normx1=norm(x,1); %求解向量 x 的 1 范数 normx1 normx2=norm(x,2); %求解向量 x 的 2 范数 normx2 normxinf=norm(x,inf); %求解向量 x 的无穷范数 normxinf normy1=norm(y,1); %求解向量 y 的 1 范数 normy1 normy2=norm(y,2); %求解向量 y 的 2 范数 normy2 normyinf=norm(y,inf); %求解向量 y 的无穷范数 normyinf z1=[normx1,normx2,normxinf]; z2=[normy1,normy2,normyinf]; end
for i=2:n
for j=i:n U(i,j)=A(i,j)-L(i,1:i-1)*U(1:i-1,j);
式
%Doolittle 分解计算上三角矩阵的公
L(j,i)=(A(j,i)-L(j,1:i-1)*U(1:i-1,i))/U(i,i); %Doolittle 分解计算下三角矩 阵的公式
end
1 1 1 ������ x = (1, 2 , 3 , … , ������) ,
������ = (1,2, … , ������)������.
对n = 10，100，1000甚至更大的n计算其范数，你会发现什么结果?你能否修改
你的程序使得计算结果相对精确呢？
1.1 源代码
function [z1,z2]=norm_vector(n) %向量 z1 的值为向量 x 的是三种范数，向量 z2 的值为向量 y 的三 种范数，n 为输入参数

matlab仿真实例100题

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实现

一、最小二乘法，用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是一款流行的科学软件，广泛应用于数学建模、数据分析、图像处理等领域。

本文将通过几个实际案例，介绍如何利用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)为截距。

  1. 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

解答：
应用原有迭代公式： 程序： clc clear x(1)=1; i=1; format while (f(x(i))~=0)
i=i+1; x(i)=sqrt(10/(x(i-1)+4)); end fprintf('方程在[1,1.5]内的解为：x=%f 迭代次数为：k=%d',x(i),i)
《矩阵与数值分析》 实验报告
第一题
题目：
一、对于数列1, 1 , 1 , 1 , 1 ,，有如下两种生成方式 3 9 27 81
1、首项为 a0
1，递推公式为 an

1 3
an1, n
1, 2, ；
2、前两项为 a0
1, a1

1 3
，递推公式为
an

10 3
an1
an2 , n
x0=x; x=B*x0+f; k=k+1; end fprintf('方程组的近似解为：') x fprintf('迭代次数为：%d 次',k) 运行结果：
Gauss-Seidel 迭代法： 程序： clc clear A=[6 2 1 -2;2 5 0 -2;-2 0 8 5;1 3 2 7]; b=[4;7;-1;0]; D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); B=inv(D-L)*U; f=inv(D-L)*b; x0=[0;0;0;0]; x=B*x0+f; k=0; while(max(abs(x(1:3)-x0(1:3)))>=1e-6)
a=abs(A(i:n,i)); d=max(a); l=find(a(1:(n-i+1))==d); l=l+i-1; c=A(l,:); A(l,:)=A(i,:); A(i,:)=c; for j=i+1:n
l1=A(j,i)/A(i,i); A(j,:)=A(j,:)-l1*A(i,:); end end x=zeros(n,1); x(n)=A(n,n+1)/A(n,n); for i=1:n-1 z=A(n-i,n+1); for j=n-i+1:n; z=z-A(n-i,j)*x(j); end x(n-i)=z/A(n-i,n-i); end fprintf(‘方程组的解为：\n’); x 运行结果：
end f=fval; subplot(2,3,k); plot(val,fval(1:length(val))); xlabel('x'); ylabel('y'); fprintf(' 在 区 间 [%f,%f] 内 \n',x(k),x(k+1)); fprintf('三次样条函数 s（%d）=',k); pretty(ans) end

ቤተ መጻሕፍቲ ባይዱ
2, 3, ；
给出利用上述两种递推公式生成的序列的第 50 项。
解答：
第一小题： 程序： clc a(1) = sym(1); for j=1:49
a(j+1)=a(j)/sym(3); end %format rat disp('第五十项为：') disp(a(50)) 运行结果：
第二小题： 程序： clc a(1)=sym(1); a(2)=sym(1/3); for j=1:50
if i~=j w=w*(x0(j)-x0(i));
end end f(k)=f(k)+y0(j)/w; end end syms x px=f(1); for i=2:n d=1; for j=2:i d=d*(x-x0(j-1)); end px=px+f(i)*d; end p=px;
调用程序： function cdd_five a=-1; b=1; n=[3 6 9 12] for i=1:4 subplot(2,2,i) x=[a:0.001:b]; plot(x,f(x)); text(0,f(0),'y=f(x)'); hold on clear y0 y0=getxy(a,b,n(i)); p=cdd_Newton(a,b,y0,n(i)); fprintf('在 n=%d 时 Newton 插值多项式为： ',n(i)); px=simplify(p); px=vpa(px); pretty(px); plot(x,subs(p,x),'r'); text(0,subs(p,0),'y=p(x)'); end
pp(k)=simplify(pp(k)); coeff=sym2poly(pp(k)); ans=sym(coeff,'d'); ans=poly2sym(ans,'t'); val=x(k):0.001:x(k+1); for i=1:length(val)
fval(i)=coeff(1)*val(i)^3+coeff(2)*val(i)^2+c oeff(3)*val(i)+coeff(4);
function y=f(x) y=x^3+4*x^2-10; 运行结果：
Aitken 加速后： 程序：
clc; k=1; x(1)=1; while(f(x(k))~=0)
k=k+1; a=sqrt(10/(x(k-1)+4));
b=sqrt(10/(a+4)); x(k)=b-(b-a)^2/(b-2*a+x(k-1)); end fprintf('方程在[1,1.5]内的根为: x=%f 迭代次数为：k=%d',x(k),k); function varout=f(x) varout=x^3+4*x^2-10; 运行结果：
0.45 0.6708
0.53 0.7280
求解，其中边界条件为
.
三次样条差值函数程序： function cdd_chazhi(x,y,,s) n=size(x,2); for i=1:n-1
h(i)=x(i+1)-x(i); end for i=2:n-1
lm(i)=h(i)/(h(i)+h(i-1)); u(i)=1-lm(i);
z=z-A(n-i,j)*x(j); end x(n-i)=z/A(n-i,n-i); end fprintf('方程组的解为：\n'); x 运行结果：
题目：
四、已知一组数据点 数 满足
第四题
，编写一程序求解三次样条插值函
并针对下面一组具体实验数据
0.25
0.3
0.5000
0.5477
0.39 0.6245
程序：
function p=cdd_Newton(a,b,y0,n) for i=1:n
x0(i)=a+(i-1)*(b-a)/(n-1); %y0(i)=fx(x0(i)); %y0(i)=1/(1+25*x0(i)^2); end for k=1:n f(k)=0; for j=1:k
w=1; for i=1:k
g(i)=3*(u(i)*(y(i+1)-y(i))/h(i)+lm(i)*(y(i)-y(i1))/y(i-1)); end g(1)=3*(y(2)-y(1))/h(1)-h(1)*d2s(1)/2; g(n)=3*(y(n)-y(n-1))/h(n-1)+h(n-1)*d2s(2)/2 lm(n)=1; u(1)=1; lm u A=diag(lm(2:n),-1)+diag(u,1); for i=1:n
题目：
[0.39,0.45]
第五题
[0.45,0.53]
五、编写程序构造区间
上的以等分结点为插值结点的 Newton 插值公式，假设结点数
为 （包括两个端点），给定相应的函数值，插值区间和等分的份数，该程序能快速计算出
相应的插值公式。以
，
为例计算其对应的插值公式，
分别取不同的 值并画出原函数的图像以及插值函数的图像，观察当 增大时的逼近效果.
function y=f(x) y=1./(1+25*x.^2); function y=getxy(a,b,n) for i=1:n
x0(i)=a+(i-1)*(b-a)/(n-1); y0(i)=f(x0(i)); end y=y0;
运行结果：
图形： 由结果可知，n 越大，插值函数曲线越逼近原函数
a(j+2)=10*a(j+1)/sym(3)-a(j); end format rat disp('第五十项为:') disp(a(50)) 运行结果：
第二题
题目：
二、利用迭代格式
xk 1
10 , k 0,1, 2, xk 4
及 Aitken 加速后的新迭代格式求方程 x3 4x2 10 0 在[1, 1.5] 内的根
迭代法计算停止的条件为：
max
1 j3
x (k 1) j

x(k) j
106 ．
2. 用 Gauss 列主元消去法、QR 方法求解如下方程组：
2

4
4
2 1 2
1 3 0
2 1 1


x1 x2 x3



1 2 1
QR 方法： 程序：
function three_QR(A) clc b=size(A); n=b(1); for i=1:n-1
Q=eye(n); a=A(i:n,i); a(1)=a(1)-sqrt(sum(a.^2)); w=a; H=eye(n-i+1)-2/(w'*w)*(w*w'); Q(i:n,i:n)=H; A=Q*A; end x=zeros(n,1); x(n)=A(n,n+1)/A(n,n); for i=1:n-1 z=A(n-i,n+1); for j=n-i+1:n;