浅谈对矩阵理论、数值分析及matlab应用的理解7

1、附录1Ax=b1、利用高斯消元法clearformatratA=randn(5)b=randn(5,1)m,n=size(A);fori=1:(m-1)numb=int2str(i);forj=(i+1):mA(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i);endAendx(m)=A(m,n)/A(m,m);fori=(m-1):-1:1x(i)=(A(i,n)-A(i,i+1:m)*x(i+1:m)/A(i,i);endx2、利用高斯列主元消元法clearformatratA=randn(5)b=randn(5,1)m,n=size(A);fori=1:(m-1)numb

2、=int2str(i);disp(第,numb,次选列主元后的增广矩阵)temp=max(abs(A(i:m,i);a,b=find(abs(A(i:m,i)=temp);tempo=A(a(1)+i-1,:);A(a(1)+i-1,:)=A(i,:);A(i,:)=tempodisp(第,numb,次消元后的增广矩阵)forj=(i+1):mA(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i);endend%回代过程disp(回代求解)x(m)=A(m,n)/A(m,m);fori=(m-1):-1:1x(i)=(A(i,n)-A(i,i+1:m)*x(i+1:m)/A(i,

3、i);endx3、LU分解A=randn(5);L,U,P=lu(A)L=PLn=5;A=randn(n);fori=1:nA(i,i)=sum(abs(A(i,:);endb=randn(n,1);D=diag(diag(A);L=tril(A,-1)U=triu(A,1)M=D(L+U);f=Db;Rho=max(abs(eig(M)R=1e-07;switchsign(1-Rho)case-1disp(TheJacobianMethodisnotapplicable)otherwisex(:,1)=normrnd(0,9,n,1);k=1;whilek=Rk=k+1;elseX=x(:,

4、k+1)IterN=kbreakendendendY=AbER_1=norm(Y-X,1)ER_2=norm(Y-X,2)ER_inf=norm(Y-X,inf)改变谱半径的大小看收敛速度的变化n=5;A=randn(n)fori=1:nA(i,i)=sum(abs(A(i,:);endb=randn(n,1)D=diag(diag(A);L=tril(A,-1);U=triu(A,1);M=D(L+U);f=Db;Rho=max(abs(eig(M)/2R=1e-07;switchsign(1-Rho)case-1disp(TheJacobianMethodisnotapplicable)o

5、therwisex(:,1)=normrnd(0,9,n,1);k=1;whilek=Rk=k+1;elseX=x(:,k+1)IterN=kbreakendendendY=Ab;ER_1=norm(Y-X,1);ER_2=norm(Y-X,2);ER_inf=norm(Y-X,inf);n=5;A=randn(n)fori=1:nA(i,i)=sum(abs(A(i,:);endb=randn(n,1)D=diag(diag(A);L=tril(A,-1);U=triu(A,1);M=D(L+U);f=Db;Rho=max(abs(eig(M)/3R=1e-07;switchsign(1-R

6、ho)case-1disp(TheJacobianMethodisnotapplicable)otherwisex(:,1)=normrnd(0,9,n,1);k=1;whilek=Rk=k+1;elseX=x(:,k+1)IterN=kbreakendendendY=Ab;ER_1=norm(Y-X,1);ER_2=norm(Y-X,2);ER_inf=norm(Y-X,inf);n=5;C=normrnd(1,4,n);AR=qr(C);A=A+3*eye(n);CondA=cond(A,2)X=ones(n,1);b=A*XY=AbD=diag(diag(A);L=tril(A,-1);

7、U=triu(A,1);fork=1:100W(k)=0.2*k;Gs(:,:,k)=(D-W(k)*L)(1-W(k)*D+W(k)*U);Rho(k)=max(abs(eig(Gs(:,:,k);endMinRho,k=min(Rho)IterN=kw=W(K)GS=(D-w*L)(1-w)*D+w*U);FS=(D-w*L)(w*b);r=1e-06;clearallx=-3:0.1:3;Fun=(x)-2*x.A6+10*x.A4-2*x.A3+5;fun=Fun(x);plot(x,fun,x,0*x,LineWidth,3)holdon2, 2.5, Respectivelytit

8、le(Solvetoalgebraicequation,fontsize,24)%RealRootsarein-1,-0.5,0.5,1anda=1;b=3;X=(a+b)/2;Beta=(sqrt(5)-1)/2;x1=a+(1-Beta)*(b-a);x2=a+Beta*(b-a);fa=Fun(a);fb=Fun(b);fx=Fun(X);f1=Fun(x1);f2=Fun(x2);xr=1e-05;fr=1e-05;N=0;whileabs(fx)=frN=N+1ifabs(b-a)xrX=(a+b)/2breakelseiff1=0X=x1breakelseiff2=0X=x2bre

9、akelseiffa*f10b=x1;X=(a+b)/2;x1=a+(1-Beta)*(b-a);x2=a+Beta*(b-a);f1=Fun(x1);f2=Fun(x2);fx=Fun(X);elseiffb*f20a=x2;X=(a+b)/2;x1=a+(1-Beta)*(b-a);x2=a+Beta*(b-a);f1=Fun(x1);f2=Fun(x2);fx=Fun(X);elsea=x1;b=x2;X=(a+b)/2;x1=a+(1-Beta)*(b-a);x2=a+Beta*(b-a);f1=Fun(x1);f2=Fun(x2);fx=Fun(X);endendplot(X,Fun

10、(X),r.,markersize,30)formatlongx0=-3- 2*x.A3 +5, x0)- 2*x.A3 +5, x0)Fsolve_x=fsolve(x)-2*x.A6+10*x.A4Fzero_x=fzero(x)-2*x.A6+10*x.A4附录4functionF=Experiment_03_Fun(x)F=x(1).A2+2*x(1)*x(2)-x(2).A2-1;x(1)A2-0.94*x(1)*x(2)-x(2)A2-9;%SolvetheEquationsbyNewtonIterationMethodF=(x,y)x.A2+b*x*y-yA2-1;x.A2+c*

11、x*y+yA2-9;Fun=F1;F2;Fx=diff(Fun,s);%求偏导Fy=diff(Fun,t);%求偏导DF=FxFy;%生成雅可比矩阵x(:,1)=2;-3;%定义初始值n=50;%最大迭代数Rx=1e-8;%偏差k=1;whileknf(:,k)=F(x(1,k),x(2,k);J=subs(DF,s,x(1,k);J=subs(J,t,x(2,k);x(:,k+1)=x(:,k)-Jf(:,k);%牛顿迭代法ifnorm(x(:,k+1)-x(:,k)RxRoot=x(:,k+1)Iter=kk=n;elsek=k+1;endend%Solvesystemofnonlinea

12、requationsS=fsolve(Experiment_03_Fun,x(:,1)R=norm(S-Root)fun=F(Root(1),Root(2)figure(2)set(gca,fontsize,24)plot(x(1,:),x(2,:),c,x(1,:),x(2,:),r.,linewidth,2,markersize,30)holdonaxis(-55-55)b=3symsst;F1=s.A2+b*s*t-tA2-1;F2=s.A2+c*s*t+tA2-9;h1=ezplot(F1)set(h1,linewidth,3,color,m)holdonh2=ezplot(F2)se

13、t(h2,linewidth,3,color,b)title(TheSystemofNonlinearEquations:x.A2+b*x*y-yA2-1=0;x.A2+c*x*y+yA2-9=0;)xlabel(fontsize24fontname黑体图3-2)附录5附程序1：clearallclc%输入插值节点XI=0,2.3,4.9,9.1,13.7,18.3,22.9,27.2;YI=22.8,22.8,22.8,20.6,13.9,11.7,11.1,11.1;%UseforPlottingCurveoftheinterpolatedfunctionx=linspace(0,29,9

14、0);%调用拉格朗日插值函数Ln=Lagrange_Fun_01(x,XI,YI);%画出插值节点plot(XI,YI,ko,LineWidth,3,markersize,10)pause(3)holdon%画出拉格朗日函数图像plot(x,Ln,-,color,1,0,0,LineWidth,3)legend(数据点，Lagrange插值曲线)title(LagrangeInterpolationPolynomialofT=T(x),fontsize,18)xlabel(湖水深，FontSize，24)ylabel(温度,FontSize,24)%简化方程式s=sym(s);f=simple

15、(Lagrange_Fun_01(s,XI,YI)y=sym2poly(f)附程序2：functionLn=Lagrange_Fun_01(x,XI,YI)ifisa(x,sym)=1;n=length(XI)-1;Ln=0;Pn=sym(ones(n+1,1);fork=1:n+1fori=1:n+1ifi=kPn(k)=Pn(k)*(x-XI(i)/(XI(k)-XI(i);elsePn(k)=Pn(k);endendLn=Ln+Pn(k)*YI(k);endelsen=length(XI)-1;L=ones(n+1,length(x);Ln=zeros(size(x);fork=1:n+

16、1fori=1:n+1ifi=kL(k,:)=L(k,:).*(x-XI(i)./(XI(k)-XI(i);elseL(k,:)=L(k,:);endendLn=Ln+L(k,:).*YI(k);endend附录6三次样条差值法clearclcclf%输入插值节点XI=0;2.3;4.9;9.1;13.7;18.3;22.9;27.2;YI=22.8;22.8;22.8;20.6;13.9;11.7;11.1;11.1;%画出插值节点plot(XI,YI,ko,linewidth,3,markersize,10)pause(3)holdon%画出三次样条插值的图象PP=csape(XI,YI,

17、1,1,0,0);%自然边界条件fnplt(PP,3,0,29)legend(数据点，三次样条插值)title(TheCubicSplineInterpolationWithEndConditions)xlabel(湖水深，FontSize,24)ylabel(温度,FontSize,24)%显示每段样条函数多项式的系数fork=1:length(XI)-1P=poly2sym(Coefs(k,:);end最小二乘法拟合clearallclcx=104121134147150177180191190204;y=100125135155170185200205210235;P,S=polyfit

18、(x,y,3);t=linspace(x(1),x(end),41);%UsethefittingpolynomialPtoevaluationPP=polyval(P,t)plot(t,PP,x,y,ro,LineWidth,3,markersize,8)set(gca,FontSize,24)legend(TheFittingCurve,TheData,4)title(CurveFittingbyLeastSquareApproximation)附录7附程序1：clearallclcsymsxt;forn=3:10;m=50;a=-0.5;b=3;Ln=Legendre_Poly(n,x)

19、;Lnt=expand(subs(Ln,(b-a)(2*t-a-b);%Lnt=expand(subs(Ln,x);T=linspace(a,b,m+1);%OptimalSquareApproximationofFunctionsFun=(1+t.A2).*cos(t/3).*exp(-(t+1)/2)+1;fun=subs(Fun,T);fori=1:n+1forj=1:n+1Phi(i,j)=int(Lnt(i).*Lnt(j),a,b);endFY(i)=int(Lnt(i)*Fun,a,b);B(i)=Phi(i,i)FY(i);endd=double(int(B*Lnt-Fun)A

20、2,-0.5,3)A0.5);ifd=0.5e-5,N=nD=dbreak;elsen=n+1;endendOptm_Poly=simple(B*Lnt)Pn=sym2poly(Optm_Poly)Optm_P=subs(Optm_Poly,T);set(gca,fontsize,24)plot(T,fun,T,Optm_P,r.-,linewidth,3)title(OptimalSquareApproximationofFunctions)legend(TheCurveofTheFunction,TheCurveofTheOSAPolynomial,1)%axissquareholdonp

21、lot(T,0*T,k,linewidth,3)附程序2：Legendre_Poly.MfunctionLn=Legendre_Poly(n,x)ifisa(x,sym)=1;ifn=0Ln=sym(1);elseLn=sym(ones(n+1,1);Ln(2)=x;fork=2:nLn(k+1)=k(2*k-1)*x*Ln(k)-k(k-1)*Ln(k-1);endendelseifn=0Ln=ones(1,length(x);elseLn=ones(n+1,length(x);Ln(2,:)=x;fork=2:nLn(k+1,:)=k(2*k-1)*x.*Ln(k,:)-k(k-1)*Ln

22、(k-1,:);endendend附录8%ByTrapezoidQuadratureclearallclcx=00.21.002.103.505.006.807.509.0011.212.0;y=-1.64-1.58-1.68-1.84-1.58-0.86-0.39-0.31-0.39-0.77-0.86;Yt=0;fork=1:10Yt=Yt+(y(1,k)+y(1,k+1)*(x(1,k+1)-x(1,k)/2;endYtplot(x,y,ko,LineWidth,3,markersize,10)pause(3)holdonh=plot(x,y);pause(5)holdonpatch(x

23、(1,1)xx(1,end),0y0,y)title(复化的梯形求积法,FontSize,24)xlim(x(1,1)-0.5,x(1,end)+0.5)xlabel(水平距离(km),FontSize,24)ylabel(深度(m),FontSize,24)%Simpsonmethodclearallclcx=00.21.002.103.505.006.807.509.0011.212.0;y=-1.64-1.58-1.68-1.84-1.58-0.86-0.39-0.31-0.39-0.77-0.86;Sn=0;fork=1:2:10Sn=Sn+(y(1,k)+4*y(1,k+1)+y(1

24、,k+2)*(x(1,k+2)-x(1,k)/6;endSnplot(x,y,ko,LineWidth,3,markersize,10)pause(3)holdonh=plot(x,y);pause(5)holdonpatch(x(1,1)xx(1,end),0y0,y)title(Simpson求积法,FontSize,24)xlim(x(1,1)-0.5,x(1,end)+0.5)xlabel(水平距离(km),FontSize,24)ylabel(深度(m),FontSize,24)MonteCarloclearallclcx=00.21.002.103.505.006.807.509.

25、0011.212.0;y=-1.64-1.58-1.68-1.84-1.58-0.86-0.39-0.31-0.39-0.77-0.86;P,S=polyfit(x,y,8);a=0;b=12.0;x=linspace(a,b,10000);y=polyval(P,x);Squad1=-quad(x)polyval(P,x),a,b)m=1e5;fork=1:150rng(shuffle)X=unifrnd(a,b,1,m);p2(1,k)=mean(polyval(P,X);endIL=-(b-a)*mean(p2(1,:)m=(1.58-1.64)/(0.2);n=(0.86-0.77)/

26、(12-11.2);PP=csape(x,y,1,1,m,n)t=linspace(x(1),x(end),1000);y=fnval(PP,t);Squad2=-quad(t)fnval(PP,t),x(1),x(end)m=1e5;fork=1:150rng(shuffle)X=unifrnd(a,b,1,m);p2(1,k)=mean(fnval(PP,X);endIP=-(b-a)*mean(p2(1,:)附录9程序zhw0%EulerMethodfortheExample:clearallclc%DefineTheleftfunctionoftheODEa=1;b=1;F=(s,x)

27、-a.*x-s.A2*sin(s).*exp(-05*s)+b;y0=-1;Tf=2*pi;Tspan=0,Tf;T,z=ode45(F,Tspan,y0);h=0.2;t=0:h:Tf;Lt=length(t);y(1)=y0;fork=2:Lty(k)=y(k-1)+h*F(t(k-1),y(k-1);endset(gca,Fontsize,24)plot(T,z,LineWidth,4,markersize,5)holdonplot(t,y,g-o,LineWidth,4,markersize,5)%4thOrderRunge_KuttaMethodForSystemsEquations

28、h=0.2;t=0:h:Tf;Y(1)=y0;fork=1:length(t)-1;K1(k)=F(t(k),Y(k);z(k)=Y(k)+2h*K1(k);K2(k)=F(t(k)+h/2,z(k);z(k)=Y(k)+2h*K2(k);K3(k)=F(t(k)+h/2,z(k);z(k)=Y(:,k)+h*K3(k);K4(k)=F(t(k)+h,z(k);Y(k+1)=Y(k)+6h*(K1(k)+2*K2(k)+2*K3(k)+K4(k);endholdonplot(t,Y,r-,LineWidth,4,markersize,5)title(EulerMethodand4R-KMeth

29、odForSystemsEquations)legend(IntegralCurve,EulerCurve,TheCurveObtainedByOrder4R-K,2)yo=z(end)ye=y(k)Y4=Y(k)clearallclc%DefineTheleftfunctionoftheODEa=1;b=1;F=(s,x)-a.*x-s.A2*sin(s).*exp(-05*s)+b;y0=unifrnd(-1,1);Tf=4*pi;Tspan=0,Tf;T,z=ode45(F,Tspan,y0);h=0.2;t=0:h:Tf;Lt=length(t);y(1)=y0;fork=2:Lty(

30、k)=y(k-1)+h*F(t(k-1),y(k-1);endsubplot(1,2,1)plot(T,z,LineWidth,4,markersize,5)holdonplot(t,y,g-o,LineWidth,4,markersize,5)%4thOrderRunge_KuttaMethodForSystemsEquationsh=0.2;t=0:h:Tf;Y(1)=y0;fork=1:length(t)-1;K1(k)=F(t(k),Y(k);z(k)=Y(k)+2h*K1(k);K2(k)=F(t(k)+h/2,z(k);z(k)=Y(k)+2h*K2(k);K3(k)=F(t(k)

31、+h/2,z(k);z(k)=Y(:,k)+h*K3(k);K4(k)=F(t(k)+h,z(k);Y(k+1)=Y(k)+6h*(K1(k)+2*K2(k)+2*K3(k)+K4(k);endholdonplot(t,Y,r-,LineWidth,4,markersize,5)title(EulerMethodand4R-KMethodForSystemsEquations)legend(IntegralCurve,EulerCurve,TheCurveObtainedByOrder4R-K,1)%EulerMethodfortheExample:Tf=10*pi;Tspan=0,Tf;T,

32、z=ode45(F,Tspan,y0);h=0.2;t=0:h:Tf;Lt=length(t);y(1)=y0;fork=2:Lty(k)=y(k-1)+h*F(t(k-1),y(k-1);endsubplot(1,2,2)plot(T,z,LineWidth,4,markersize,5)holdonplot(t,y,g-o,LineWidth,4,markersize,5)%4thOrderRunge_KuttaMethodForSystemsEquationsh=0.2;t=0:h:Tf;Y(1)=y0;fork=1:length(t)-1;K1(k)=F(t(k),Y(k);z(k)=

33、Y(k)+2h*K1(k);K2(k)=F(t(k)+h/2,z(k);z(k)=Y(k)+2h*K2(k);K3(k)=F(t(k)+h/2,z(k);z(k)=Y(:,k)+h*K3(k);K4(k)=F(t(k)+h,z(k);Y(k+1)=Y(k)+6h*(K1(k)+2*K2(k)+2*K3(k)+K4(k);endholdonplot(t,Y,r-,LineWidth,4,markersize,5)title(EulerMethodand4R-KMethodForSystemsEquations)legend(IntegralCurve,EulerCurve,TheCurveObt

34、ainedByOrder4R-K,1)附录10解clearallclcTf=16*pi;y0=0;0;subplot(3,2,1)a=2;b=0.75;c=1;d=-0.5;Tspan=0,Tf;Options=odeset(RelTol,1e-5,AbsTol,1e-6);T,Z=ODE_Parameter_Fun(a,b,c,d,Tf,y0)plot(T,Z,LineWidth,2.5)holdonplot(T,0*T,k,LineWidth,2.5)subplot(3,2,2)a=0.1;b=9.0025;c=0;S=roots(1;a;b);%求方程的根d=S(1,1)%取第一个根Ts

35、pan=0,Tf;Options=odeset(RelTol,1e-5,AbsTol,1e-6);T,Z=ODE_Parameter_Fun(a,b,c,d,Tf,y0)plot(T,Z,LineWidth,2.5)holdonplot(T,0*T,k,LineWidth,2.5)subplot(3,2,3)a=0;b=9;c=0;S=roots(1;a;b);d=0.95*S(1,1)Tspan=0,Tf;Options=odeset(RelTol,1e-5,AbsTol,1e-6);T,Z=ODE_Parameter_Fun(a,b,c,d,Tf,y0)plot(T,Z,LineWidth

36、,2.5)holdonplot(T,0*T,k,LineWidth,2.5)subplot(3,2,4)a=0.1;b=9.0025;c=0;S=roots(1;a;b);d=S(1,1)Tspan=0,Tf;Options=odeset(RelTol,1e-5,AbsTol,1e-6);T,Z=ODE_Parameter_Fun(a,b,c,d,Tf,y0)plot(T,Z,LineWidth,2.5)holdonplot(T,0*T,k,LineWidth,2.5)subplot(3,2,5)a=2;b=1.25;c=1;d=0;Tspan=0,Tf;Options=odeset(RelT

37、ol,1e-5,AbsTol,1e-6);T,Z=ODE_Parameter_Fun(a,b,c,d,Tf,y0)plot(T,Z,LineWidth,2.5)holdonplot(T,0*T,k,LineWidth,2.5)functionT,Z=ODE_Parameter_Fun5(a,b,c,d,Tf,y0)Tspan=0,Tf;T,Z=ode45(ODE_Parameter,Tspan,y0);functiondydt=ODE_Parameter(t,y)y=y(1);y(2);dydt=y(2)-a*y(2)-b*y(1)+c+cos(imag(d)*t)*exp(real(d)*t

38、);endend解clearallclcTf=20*pi;y0=0;0;symsy(t)Dy=diff(y,1);a=2;b=0.75;c=1;d=-0.5;Y1=dsolve(diff(y,2)=-a*diff(y,1)-b*y+c+cos(imag(d)*t)*exp(real(d)*t),y(0)=0,Dy(0)=0)a=0.1;b=9.0025;c=0;S=roots(1;a;b);d=S(1,1);Y2=dsolve(diff(y,2)=-a*diff(y,1)-b*y+c+cos(imag(d)*t)*exp(real(d)*t),y(0)=0,Dy(0)=0)a=0;b=9;c=

39、0;S=roots(1;a;b);d=0.95*S(1,1);Y3=dsolve(diff(y,2)=-a*diff(y,1)-b*y+c+cos(imag(d)*t)*exp(real(d)*t),y(0)=0,Dy(0)=0)a=0.1;b=9.0025;c=0;S=roots(1;a;b);d=S(1,1);Y4=dsolve(diff(y,2)=-a*diff(y,1)-b*y+c+cos(imag(d)*t)*exp(real(d)*t),y(0)=0,Dy(0)=0)a=2;b=1.25;c=1;d=0;Y5=dsolve(diff(y,2)=-a*diff(y,1)-b*y+c+

40、cos(imag(d)*t)*exp(real(d)*t),y(0)=0,Dy(0)=0)附录11vdp:functionxprime=vdp(t,x)globalmu;xprime=x(2);-x(1)+mu*(1-x(12)*x(2);clearclcy0=-3,0.3;Tf=80;Tspan=0,Tf;globalmu;m=-1-0.8-0.6-0.4-0.200.20.40.60.81;i=1;whilei=11mu=m(i);t,x=ode45(vdp,Tspan,y0);figure(1)subplot(3,4,i)plot(t,x,t,0*t,k,LineWidth,2)xlabel(Time)titleName=strcat(mu=,num2str(mu);title(titleName);figure(2)subplot(3,4,i