MATLAB数学建模与仿真（第2版·微课视频版）程序代码 第7章 MATLAB科学计算

【例7-1】使用除法求解系数矩阵可逆的恰定线性方程组。在命令窗口中输入如下语句：A=pascal(4)det_A=det(A)B=rand(4,1)X1=A\BX2=inv(A)*B命令窗口中的输出结果如下所示：A=1111123413610141020det_A=1B=0.81470.90580.12700.9134X1=-2.58139.1360-8.17512.4351X2=-2.58139.1360-8.17512.4351【例7-2】使用除法求解系数矩阵不可逆的恰定线性方程组。在命令窗口中输入如下语句：A=[137;-144;11018];det_A=det(A)B=[6;4;15]X=A\B命令窗口中的输出结果如下所示：det_A=0B=6415Warning:Matrixissingulartoworkingprecision.X=NaNInf-Inf【例7-3】使用除法求解欠定线性方程组。在命令窗口中输入如下语句：C=magic(4);A=C(1:3,:)B=[1;0;0];X=A\B命令窗口中的输出结果如下所示：A=16231351110897612X=0.18630.02940-0.1569【例7-4】使用除法求解超定线性方程组。在命令窗口中输入如下语句：T=magic(5)A=T(:,1:4)B=[1;0;0;0;0];X=A\B命令窗口中的输出结果如下所示：T=17241815235714164613202210121921311182529A=172418235714461320101219211118252X=-0.00410.0437-0.03050.0060【例7-5】使用伪逆矩阵的方法求解奇异矩阵的线性方程的解。在命令窗口中输入如下语句：A=[137;-144;11018];B=[5;2;12];X=pinv(A)*BC=A*X命令窗口中的输出结果如下所示：X=0.3850-0.11030.7066C=5.00002.000012.0000【例7-6】使用求逆法计算线性方程组的所有解。在命令窗口中输入如下语句：A=[1234;5678;9101112];B=[1;1;2];X1=null(A)X2=pinv(A)*B命令窗口中的输出结果如下所示：X1=0.51350.1906-0.82670.12870.1129-0.82900.20030.5098X2=-0.1250-0.02080.08330.1875【例7-7】使用LU分解法求解线性方程组。在命令窗口中输入如下语句：A=[211;1-23;6-51];B=[1;5;7];C=det(A)[L,U]=lu(A)X=U\(L\B)命令窗口中的输出结果如下所示：C=50L=0.33331.000000.1667-0.43751.00001.000000U=6.0000-5.00001.000002.66670.6667003.1250X=0.3600-0.76001.0400【例7-8】使用LU分解法求解线性方程组。在命令窗口中输入如下语句：A=[-185;9-12;2-57];B=[2;3;5];C=det(A)[L,U]=lu(A)X=U\(L\B)命令窗口中的输出结果如下所示：C=-690L=-0.11111.000001.0000000.2222-0.60561.0000U=9.0000-1.00002.000007.88895.2222009.7183X=0.1913-0.09570.5913【例7-9】使用Cholesky分解法求解线性方程组。在命令窗口中输入如下语句：A=[213;154;346];B=[1;3;5];C=det(A)R=chol(A)TR=R'X=R\(TR\B)命令窗口中的输出结果如下所示：C=1R=1.41420.70712.121302.12131.1785000.3333TR=1.4142000.70712.121302.12131.17850.3333X=-23.0000-10.000019.0000【例7-10】对对称正定矩阵进行Cholesky分解。在命令窗口中输入如下语句：n=5;X=pascal(n)R=chol(X)C=transpose(R)*R命令窗口中的输出结果如下所示：X=111111234513610151410203515153570R=1111101234001360001400001C=111111234513610151410203515153570由结果可以看出，R是上三角矩阵，同时满足，表明上面的Cholesky分解过程是正确的。如果修改矩阵的信息，在命令窗口中继续输入如下语句：X(n,n)=X(n,n)-1[R1,p]=chol(X)C1=transpose(R1)*R1C2=X(1:p-1,1:p-1)命令窗口中的输出结果如下所示：X=111111234513610151410203515153569R1=1111012300130001p=5C1=1111123413610141020C2=1111123413610141020【例7-11】使用QR分解法求解线性方程组。在命令窗口中输入如下语句：A=[2134;1542;346-5];B=[1;3;5];[Q,R]=qr(A)X=R\(Q\B)命令窗口中的输出结果如下所示：Q=-0.53450.4257-0.7301-0.2673-0.9047-0.3319-0.80180.01770.5974R=-3.7417-5.0780-7.48331.33630-4.0267-2.2351-0.1951000.0664-6.5709X=00.28460.4878-0.1870【例7-12】使用共轭梯度法求解线性方程组。在命令窗口中输入如下语句：A=[311;2-15;8-40];B=[2;4;1];[X,flag,relres,iter,resvec]=bicg(A,B)命令窗口中的输出结果如下所示：X=0.30000.35000.7500flag=0relres=5.4820e-016iter=3resvec=4.58263.93706.60520.0000其中，flag表示在默认迭代次数内收敛，relres表示相对残差norm(B-A*x)/norm(B)，iter表示终止的迭代次数，resvec表示每次迭代的残差。【例7-13】计算一元函数在[-3,4]区间的零点。首先绘制函数的曲线，在命令窗口中输入如下语句：x=-3:0.1:4;y=x.*x.*sin(x)-x+1;plot(x,y,'r')xlabel('x');ylabel('f(x)');title('Thezerooffunction')holdonh=line([-3,4],[0,0]);set(h,'color','g')grid;图形窗口中的输出结果如图7-1所示。图7-1一元函数曲线在求解函数零点之前，需要绘制函数的图形，是为了在后面的步骤中使用fzero命令时，更加方便地选择初始数值0。由图7-1不难看出，曲线在[-3,4]区间内包含3个零点。其次计算函数某点附近的零点，在命令窗口中输入如下语句：f=@(x)x-10*sin(x);X1=fzero(f,-2.5)X2=fzero(f,-1.5)X3=fzero(f,3)命令窗口中的输出结果如下所示：X1=-2.8523X2=-2.8523X3=2.8523【例7-14】求解二元方程组的零点。首先绘制函数的曲线，在命令窗口中输入如下语句：x=[-5:0.1:5];y=x;[X,Y]=meshgrid(x,y);Z=2*X-Y-exp(-X);surf(X,Y,Z)xlabel('x')ylabel('y')zlabel('z')title('Thefigureofthefunction')图形窗口中的输出结果如图7-2所示。图7-2二元方程组曲线编写求解函数的M文件，在其中输入如下的程序代码：functionF=myfun8_15(x)F=[2*x(1)-x(2)-exp(-x(1));-x(1)+2*x(2)-exp(-x(2))];将上述程序代码保存为myfun8_15.m文件。下面求解二元函数的零点，在命令窗口中输入如下语句：x0=[-5;-5];options=optimset('Display','iter');x=fsolve(@myfun8_15,x0,options)命令窗口中的输出结果如下所示：NormofFirst-orderTrust-regionIterationFunc-countf(x)stepoptimalityradius0347071.22.29e+0411612003.415.75e+031293147.0211.47e+031312854.45213881415239.5271107151867.0412130.8162116.704219.0517242.4278812.2618270.0326580.7595110.2062.59307.03149e-060.1119270.002942.510333.29525e-130.001691326.36e-072.5Equationsolved.fsolvecompletedbecausethevectoroffunctionvaluesisnearzeroasmeasuredbythevalueofthefunctiontolerance,andtheproblemappearsregularasmeasuredbythegradient.<stoppingcriteriadetails>x=0.56710.5671【例7-15】求解。首先对非线性方程组进行函数描述，并保存为myfun7_15.m，其内容如下所示：functionT=myfun7_15(x)T=x^4-[1,7;-11,9];其次对非线性方程组进行求解，在命令窗口中输入如下语句：x0=[31;21];options=optimset('Display','off');X=fsolve(@myfun8_15,x0,options)命令窗口中的输出结果如下所示：X=1.46930.3875-0.60891.9121【例7-16】已知微分方程,该方程为显式常微分方程，分别取和求解该方程。首先对微分方程进行变换得到如下形式：其次对方程组进行函数描述，并保存为myfun7_16.m，其内容如下所示：functionoutput=myfun7_16(t,y,mu)output=zeros(2,1);output(1)=y(2);output(2)=mu*(1-y(1)^2)*y(2)-y(1);再次对方程组进行求解，在命令窗口中输入如下语句：[t1,y1]=ode45(@myfun7_16,[030],[0;2],[],3);%mu=3[t2,y2]=ode45(@myfun7_16,[030],[0;2],[],5);%mu=5plot(t1,y1(:,1),'-',t2,y2(:,2),'--')title('显式常微分方程的解');xlabel('t');ylabel('y');legend('mu=3','mu=5');图形窗口中的输出结果如图7-3所示。图7-3显式常微分方程的解【例7-17】求解微分方程，该方程为线性隐式常微分方程。首先根据微分方程和通式，得到：其次对进行函数描述，并保存为myfun7_17f.m，其内容如下所示：functionoutput=myfun7_17f(t,y)output=3*y.^3+y+4;再次对进行函数描述，并保存为myfun8_17M.m，其内容如下所示：functionoutput=myfun7_17M(t,y)output=t*y.^2+1;最后对方程进行求解，在命令窗口中输入如下语句：options=odeset('RelTol',1e-6,'OutputFcn','odeplot','Mass',@myfun7_17M);[t,y]=ode45(@myfun7_17f,[010],2,options);xlabel('t');ylabel('y');title('线性隐式常微分方程的解')图形窗口中的输出结果如图7-4所示。图7-4线性隐式常微分方程的解【例7-18】求解微分方程该方程为完全隐式常微分方程。首先对方程进行函数描述，并保存为myfun7_18.m，其内容如下所示：functionoutput=myfun7_18(t,y,dydt)output=t*y.^2*dydt.^3-2*y.^3*dydt.^2+3*t*(t^2+1)*dydt-t^2*y;其次对方程进行求解，在命令窗口中输入如下语句：t0=1;y0=sqrt(3/2);yp0=0;[y0,yp0]=decic(@myfun7_18,t0,y0,1,yp0,0);[t,y]=ode15i(@myfun7_18,[120],y0,yp0);plot(t,y);xlabel('t');ylabel('y');title('完全隐式常微分方程的解');图形窗口中的输出结果如图7-5所示。图7-5完全隐式常微分方程的解【例7-19】生成[15]上均匀分布的的随机数矩阵。在命令窗口中输入如下语句：x=unifrnd(1,5,4,4)命令窗口中的输出结果如下所示：x=3.52944.83004.82872.68701.39024.85962.94154.66292.11401.63054.20114.16883.18754.88241.56754.8380【例7-20】应用求最大值、最小值的函数。在命令窗口中输入如下语句：x=1:20;y=randn(1,20);figure;holdon;plot(x,y);[y_max,I_max]=max(y)%求向量最大值及其对应下标plot(x(I_max),y_max,'ro');[y_min,I_min]=min(y)%求向量最小值及其对应下标plot(x(I_min),y_min,'g*');xlabel('x');ylabel('y');legend(‘原始数据’,‘最大值’,‘最小值’);命令窗口中的输出结果如下所示：y_max=1.5326I_max=16y_min=-1.2141I_min=13图形窗口中的输出结果如图8-6所示。图7-6最大值与最小值【例7-21】应用求和与求积的函数。在命令窗口中输入如下语句：x=1:20;y=randn(1,20);y_sum=sum(y)y_prod=prod(y)命令窗口中的输出结果如下所示：y_sum=14.298498925524077y_prod=-0.009046699652797【例7-22】应用求平均值和中值的函数，在命令窗口中输入如下语句：x=1:20;y=randn(1,20);y_mean=mean(y)%求向量平均值y_median=median(y)%求向量中间值命令窗口中的输出结果如下所示：y_mean=0.0283y_median=0.1461【例7-23】计算向量x的标准差。在命令窗口中输入如下语句：x=1:20;mean_x=mean(x);r=0;fori=1:20r=r+(x(i)-mean_x)^2;endr1=sqrt(r/20)r2=sqrt(r/19)r3=std(x)命令窗口中的输出结果如下所示：r1=5.766281297335398r2=5.916079783099616r3=5.916079783099616【例7-24】计算协方差与相关系数。在命令窗口中输入如下语句：X=[1013]';Y=[-3515]';A1=cov(X)A2=cov(X,Y)A3=corrcoef(X)A4=corrcoef(X,Y)命令窗口中的输出结果如下所示：A1=1.5833A2=1.58331.00001.000014.6667A3=1A4=1.00000.20750.20751.0000【例7-25】对随机产生的矩阵A分别进行降序和升序排序。在命令窗口中输入如下语句：A=unifrnd(1,2,4,5)Y1=sort(A,'descend')Y2=sort(A,'ascend')命令窗口中的输出结果如下所示：A=1.79431.60201.74821.91331.99611.31121.26301.45051.15241.07821.52851.65411.08381.82581.44271.16561.68921.22901.53831.1067Y1=1.79431.68921.74821.91331.99611.52851.65411.45051.82581.44271.31121.60201.22901.53831.10671.16561.26301.08381.15241.0782Y2=1.16561.26301.08381.15241.07821.31121.60201.22901.53831.10671.52851.65411.45051.82581.44271.79431.68921.74821.91331.9961【例7-26】利用一维快速傅立叶插值实现数据增采样，代码设置如下：%interpft_example.m%一维快速傅立叶插值实现数据增采样x=0:1.2:10;y=sin(x);n=2*length(x);%增采样1倍yi=interpft(y,n);%一维快速傅立叶插值xi=0:0.6:10.4;holdon;plot(x,y,'ro');%画图plot(xi,yi,'b.-');title('一维快速傅立叶插值');legend('原始数据','插值结果');由上述程序可生成图7-9所示结果。图7-9一维快速傅立叶插值新建例子新建例子【例7-27】求归一化高斯函数在区间[-11]上的定积分，代码设置如下：%quad_exam.m%求归一化高斯函数在区间[-11]上的定积分，并求得到积分过程的中间节点y=@(x)1/sqrt(pi)*exp(-x.^2);%归一化高斯函数integral(y,-1,1)%求定积分，并显示中间迭代过程fplot(y,[-11],'b');%画出函数holdon;%跟踪数据(运行完上面程序后，可以在命令行复制这些数据)trace=[9-1.00000000005.43160000e-0010.1804679399;11-1.00000000002.71580000e-0010.0728222057;13-0.72842000002.71580000e-0010.1076454255;15-0.45684000009.13680000e-0010.4817487615;17-0.45684000004.56840000e-0010.2408826755;19-0.45684000002.28420000e-0010.1142172651;21-0.22842000002.28420000e-0010.1266655031;230.00000000004.56840000e-0010.2408826755;250.00000000002.28420000e-0010.1266655031;270.22842000002.28420000e-0010.1142172651;290.45684000005.43160000e-0010.1804679399;310.45684000002.71580000e-0010.1076454255;330.72842000002.71580000e-0010.0728222057];x1=trace(:,2);%积分过程的中间节点y1=y(x1);%中间节点的函数值plot(x1,y1,'ro');%画图xlabel('x');ylabel('y');legend('高斯函数','求积分过程的中间节点');由上述代码得到输出结果如下：ans=0.8427积分过程如图7-14所示。图7-14积分过程函数quadl()的调用格式与函数quad()相同。在一维积分的过程中，有可能出现三种警告信息：'Minimumstepsizereached'，意味着子区间的长度与计算机舍入误差相当，无法继续计算了。原因可能是有不可积的奇点'Maximumfunctioncountexceeded'，意味着积分递归计算超过了10000次。原因可能是有不可积的奇点'InfiniteorNot-a-Numberfunctionvalueencountered'，意味着在积分计算时，区间内出现了浮点数溢出或者被零除。新增例子新增例子【例7-28】求，就可以用矢量积分，具体代码如下：%quadv_exam.m%求矢量积分y=@(x,n)1./(sqrt(2*pi).*(1:n)).*exp(-x.^2./(2*(1:n).^2));%归一化高斯函数integral(@(x)y(x,5),-1,1,'ArrayValued',true)由上述语句得到的输出结果如下：ans=0.68270.38290.26110.19740.1585矢量积分的结果是一个向量，其每一元素的值为一个一元函数定积分的值。新增例子新增例子【例7-29】计算二维高斯函数在矩形区间[-11-11]上的二重积分，代码设置如下：%dblquad_exam.m%计算二维高斯函数在矩形区间[-11-11]上的二重积分f=@(x,y)1/sqrt(pi)*exp(-x.^2)*1/sqrt(pi)*exp(-y.^2);%归一化高斯函数integral2(f,-1,1,-1,1)由上述语句得到输出结果如下：ans=0.7101函数integral2()处理的都是矩形积分区域。若要计算非矩形的积分区间的二重积分，可以用一个大的矩形积分区域包含非矩形的积分区间，然后在非矩形的积分区间之外的区域上把二元函数的值取零。新增例子新增例子【例7-30】计算二维高斯函数在圆形区域上的二重积分，代码设置如下：%dblquad_exam2.m%计算二维高斯函数在圆形区域sqrt(x.^2+y.^2)<1上的二重积分%在圆形区域外填充0的归一化高斯函数f=@(x,y)(1/sqrt(pi)*exp(-x.^2)*1/sqrt(pi)*exp(-y.^2)).*(sqrt(x.^2+y.^2)<=1);integral2(f,-2,2,-2,2)由上述语句得到的输出结果如下：ans=0.6321【例7-31】求解正弦函数在x0=3附近的极小值点。在命令窗口中输入如下语句：clearclcX=fminsearch(@sin,3)在命令窗口中的输出结果如下：X=4.1724【例7-32】求解约束条件下函数的极小值点。clearclcX=fmincon(@(x)3*sin(x(1))+exp(x(2)),[1;1],[],[],[],[],[00])命令窗口中的输出结果如下：Localminimumfoundthatsatisfiestheconstraints.Optimizationcompletedbecausetheobjectivefunctionisnon-decreasinginfeasibledirections,towithinthevalueoftheoptimalitytolerance,andconstraintsaresatisfiedtowithinthevalueoftheconstrainttolerance.<stoppingcriteriadetails>X=1.0e-05*0.06670.2000【例7-28】计算的线性规划。在命令窗口中输入如下语句：H=[1-1;-12];f=[-3;-5];A=[13;-12];b=[2;8];[x,fval]=quadprog(H,f,A,b,[],[],[0;0])命令窗口中的输出结果如下：Minimumfoundthatsatisfiestheconstraints.Optimizationcompletedbecausetheobjectivefunctionisnon-decreasinginfeasibledirections,towithinthevalueoftheoptimalitytolerance,andconstraintsaresatisfiedtowithinthevalueoftheconstrainttolerance.<stoppingcriteriadetails>x=1.29410.2353fval=-4.4706【例7-29】计算的线性规划在MATLAB的命令窗口中输入如下语句：f=[1;3;-5];A=[1-11;301];b=[8;4];x=linprog(f,A,b,[],[],[0;0;0],[inf;inf;inf])则命令窗口中显示输出结果如下：Optimalsolutionfound.x=004【例7-30】通过问题比较lsqnonneg函数解与无约束最小二乘解的不同。在命令窗口中输入如下语句：C=[0.03720.2869;0.68610.7071;0.62330.6245;0.63440.6170];d=[0.8587;0.1781;0.0747;0.8405];x1=lsqnonneg(C,d)x2=C\d则命令窗口显示输出结果如下：x1=00.6929x2=-2.56273.1108【例7-31】求解下列问题的最小二乘解。在命令窗口中输入如下代码：C=[0.95010.76200.61530.4057;0.23110.45640.79190.9354;0.60680.01850.92180.9169;0.48590.82140.73820.4102;0.89120.44470.17620.8936];d=[0.0578;0.3528;0.8131;0.0098;0.1388];A=[0.20270.27210.74670.4659;0.19870.19880.44500.4186;0.60370.01520.93180.8462];b=[0.5251;0.2026;0.6721];lb=-0.1*ones(4,1);ub=2*ones(4,1);x=lsqlin(C,d,A,b)命令窗口中的输出结果如下：Minimumfoundthatsatisfiestheconstraints.Optimizationcompletedbecausetheobjectivefunctionisnon-decreasinginfeasibledirections,towithinthevalueoftheoptimalitytolerance,andconstraintsaresatisfiedtowithinthevalueoftheconstrainttolerance.<stoppingcriteriadetails>x=0.1299-0.57570.42510.2438【例7-32】求解x，使得下式最小化。，初值为x=[0.30.4]由于lsqnonlin函数假设用户提供的平方和不是显式表达的，因此，lsqnonlin函数中的函数应变换为向量值函数，如下所示：编辑M文件并保存为myfun1.m文件，程序代码如下所示：functionF=myfun1(x)k=1:10;F=2+2*k-exp(k*x(1))-exp(k*x(2));其次在命令窗口中输入如下语句：x0=[0.30.4][x,resnorm]=lsqnonlin(@myfun1,x0)则命令窗口中的输出结果如下：x0=0.30000.4000Localminimumpossible.lsqnonlinstoppedbecausethesizeofthecurrentstepislessthanthevalueofthestepsizetolerance.<stoppingcriteriadetails>x=0.25780.2578resnorm=124.3622【例7-33】某玩具厂制作两种不同的涂料产片A和B，已知生产A涂料100kg需要8个工时，生产B涂料100kg需要10个工时。限定每日的工时数为40，并希望不需要临时工，也不需要工人加班生产。这两种涂料每100kg可获利100元。此外，有个顾客需求供应B涂料600kg，请问应如何制定生产计划，达到最优？分析：假设制作A和B两种涂料的数量分别为x1、x2(均以100kg计)，为了使生产计划比较合理并用人尽量少，使利润最大化，并且B涂料的产量尽量多，由以上分析可建立如下所示的数字描述：编写目标函数的M文件，保存goal.m，返回目标计算值：functio