matlab作业

实验1-1 MATLAB 运算基础1. 先求下列表达式的值，然后显示MATLAB 工作空间的使用情况并保存全部变量。

(1) 0122sin851z e =+(2) 221ln(1)2z x x =++，其中2120.455i x +⎡⎤=⎢⎥-⎣⎦(3) 0.30.330.3sin(0.3)ln , 3.0, 2.9,,2.9,3.022a a e e az a a --+=++=--(4) 2242011122123t t z t t t t t ⎧≤<⎪=-≤<⎨⎪-+≤<⎩，其中t =0:0.25:2.5 2. 已知：求下列表达式的值：(1) A+6*B 和A-B+I （其中I 为单位矩阵） (2) A*B 和A.*B (3) A^3和A.^3 (4) A/B 及B\A(5) [A,B]和[A([1,3],:);B^2]3. 设有矩阵A 和B123453166789101769,111213141502341617181920970212223242541311A B ⎡⎤⎡⎤⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥==-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦(1) 求它们的乘积C 。

(2) 将矩阵C 的右下角3×2子矩阵赋给D 。

(3) 查看MATLAB 工作空间的使用情况。

1234413134787,2033657327A B --⎡⎤⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦4. 完成下列操作：(1) 求[100,999]之间能被21整除的数的个数。

(2) 建立一个字符串向量，删除其中的大写字母。

5. 已知：⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=987654321a ，分别计算a 的数组平方和矩阵平方，并观察其结果。

6. ⎥⎦⎤⎢⎣⎡-=463521a ，⎥⎦⎤⎢⎣⎡-=263478b ，观察a 与b 之间的六种关系运算的结果。

7. 角度[]604530=x ，求x 的正弦、余弦、正切和余切。

有限元大作业程序设计学校：天津大学院系：建筑工程与力学学院专业：01级工程力学姓名：刘秀学号：\\\\\\\\\\\指导老师：连续体平面问题的有限元程序分析[题目]：如图所示的正方形薄板四周受均匀载荷的作用，该结构在边界上受正向分布压力，m kNp 1=，同时在沿对角线y 轴上受一对集中压力，载荷为2KN ，若取板厚1=t ，泊松比0=v 。

[分析过程]：由于连续平板的对称性，只需要取其在第一象限的四分之一部分参加分析，然后人为作出一些辅助线将平板“分割”成若干部分，再为每个部分选择分析单元。

采用将此模型化分为4个全等的直角三角型单元。

利用其对称性，四分之一部分的边界约束，载荷可等效如图所示。

[程序原理及实现]：用FORTRAN程序的实现。

由节点信息文件NODE.IN和单元信息文件ELEMENT.IN，经过计算分析后输出一个一般性的文件DATA.OUT。

模型基本信息由文件为BASIC.IN生成。

该程序的特点如下：问题类型：可用于计算弹性力学平面问题和平面应变问题单元类型：采用常应变三角形单元位移模式：用用线性位移模式载荷类型：节点载荷，非节点载荷应先换算为等效节点载荷材料性质：弹性体由单一的均匀材料组成约束方式：为“0”位移固定约束，为保证无刚体位移，弹性体至少应有对三个自由度的独立约束方程求解：针对半带宽刚度方程的Gauss消元法输入文件：由手工生成节点信息文件NODE.IN，和单元信息文件ELEMENT.IN结果文件：输出一般的结果文件DATA.OUT程序的原理如框图：（1）主要变量：ID：问题类型码，ID＝1时为平面应力问题，ID=2时为平面应变问题N_NODE：节点个数N_LOAD：节点载荷个数N_DOF：自由度，N_DOF=N_NODE*2（平面问题）N_ELE：单元个数N_BAND：矩阵半带宽N_BC：有约束的节点个数PE：弹性模量PR：泊松比PT：厚度LJK_ELE(I,3)：单元节点编号数组，LJK_ELE(I,1)，LJK_ELE(I,2)，LJK_ELE(I,3)分别放单元I的三个节点的整体编号X(N_NODE), Y(N_NODE)：节点坐标数组，X(I),Y(I)分别存放节点I的x，y 坐标值P_LJK(N_BC,3)：节点载荷数组，P_LJK(I,1)表示第I个作用有节点载荷的节点的编号，P_LJK(I,2)，P_LJK(I,3)分别为该节点沿x,y方向的节点载荷数值AK(N_DOF,N_BAND)：整体刚度矩阵AKE(6,6)：单元刚度矩阵BB(3,6)：位移……应变转换矩阵（三节点单元的几何矩阵）DD(3,3)：弹性矩阵SS(3,6)；应力矩阵RESULT_N(N_NOF)：节点载荷数组，存放节点载荷向量，解方程后该矩阵存放节点位移DISP_E(6):：单元的节点位移向量STS_ELE(N_ELE,3)：单元的应力分量STS_ND(N_NODE,3)：节点的应力分量（2）子程序说明：READ_IN：读入数据BAND_K：形成半带宽的整体刚度矩阵FORM_KE：计算单元刚度矩阵FORM_P：计算节点载荷CAL_AREA：计算单元面积DO_BC：处理边界条件CLA_DD：计算单元弹性矩阵SOLVE：计算节点位移CLA_BB：计算单元位移……应变关系矩阵CAL_STS：计算单元和节点应力（3）文件管理：源程序文件：chengxu.for程序需读入的数据文件：BASIC.IN，NODE.IN，ELEMENT.IN（需要手工生成）程序输出的数据文件：DATA.OUT（4）数据文件格式：需读入的模型基本信息文件BASIC.IN的格式如下表需读入的节点信息文件NODE.IN的格式如下表需读入的单元信息文件ELEMENT.IN的格式如下表输出结果文件DATA.OUT格式如下表[算例原始数据和程序分析]：（1）模型基本信息文件BASIC.IN的数据为1,4,6,5,31.,0.,1.1,1,0,2,1,0,4,1,1,5,0,1,6,0,11,-0.5,-1.5,3.,-1.,-1,6,-0.5,-0.5（2）手工准备的节点信息文件NODE.IN的数据为1 0.0 2.02 0.0 1.03 1.0 1.04 0. 0.5 1.0 0.6 2.0 0.（3）手工准备的单元信息文件ELEMENT.IN的数据为1 2 3 3 0 0 0 0 1 1 1 1 0 12 4 5 5 0 0 0 0 1 1 1 1 0 25 3 2 2 0 0 0 0 1 1 1 1 0 33 5 6 6 0 0 0 0 1 1 1 1 04 （4）源程序文件chengxu.for为：PROGRAM FEM2DDIMENSION IJK_ELE(500,3),X(500),Y(500),IJK_U(50,3),P_IJK(50,3),&RESULT_N(500),AK(500,100)D IMENSION STS_ELE(500,3),STS_ND(500,3)OPEN(4,FILE='BASIC.IN')OPEN(5,FILE='NODE.IN')OPEN(6,FILE='ELEMENT.IN')OPEN(8,FILE='DATA.OUT')OPEN(9,FILE='FOR_POST.DAT')READ(4,*)ID,N_ELE,N_NODE,N_BC,N_LOADIF(ID.EQ.1)WRITE(8,20)IF(ID.EQ.2)WRITE(8,25)20 FORMAT(/5X,'=========PLANE STRESS PROBLEM========')25 FORMAT(/5X,'=========PLANE STRAIN PROBLEM========')CALL READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,PT, & IJK_ELE,X,Y,IJK_U,P_IJK)CALL BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,& IJK_ELE,X,Y,PE,PR,PT,AK)CALL FORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK, & RESULT_N)CALL DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N)CALL SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N)CALL CAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N, & STS_ELE,STS_ND)c to putout a data fileWRITE(9,70)REAL(N_NODE),REAL(N_ELE)70 FORMAT(2f9.4)WRITE(9,71)(X(I),Y(I),RESULT_N(2*I-1),RESULT_N(2*I),& STS_ND(I,1),STS_ND(I,2),STS_ND(I,3),I=1,N_NODE)71 FORMA T(7F9.4)WRITE(9,72)(REAL(IJK_ELE(I,1)),REAL(IJK_ELE(I,2)),&REAL(IJK_ELE(I,3)),REAL(IJK_ELE(I,3)),&STS_ELE(I,1),STS_ELE(I,2),STS_ELE(I,3),I=1, N_ELE)72 FORMAT(7f9.4)cCLOSE(4)CLOSE(5)CLOSE(6)CLOSE(8)CLOSE(9)E NDcc to get the original data in order to model the problemSUBROUTINE READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR, &PT,IJK_ELE,X,Y,IJK_U,P_IJK)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),IJK_U(N_BC,3), & P_IJK(N_LOAD,3),NE_ANSYS(N_ELE,14)REAL ND_ANSYS(N_NODE,3)READ(4,*)PE,PR,PTREAD(4,*)((IJK_U(I,J),J=1,3),I=1,N_BC)READ(4,*)((P_IJK(I,J),J=1,3),I=1,N_LOAD)READ(5,*)((ND_ANSYS(I,J),J=1,3),I=1,N_NODE)READ(6,*)((NE_ANSYS(I,J),J=1,14),I=1,N_ELE)DO 10 I=1,N_NODEX(I)=ND_ANSYS(I,2)Y(I)=ND_ANSYS(I,3)10 CONTINUEDO 11 I=1,N_ELEDO 11 J=1,3IJK_ELE(I,J)=NE_ANSYS(I,J)11 CONTINUEN_BAND=0DO 20 IE=1,N_ELEDO 20 I=1,3DO 20 J=1,3IW=IABS(IJK_ELE(IE,I)-IJK_ELE(IE,J))IF(N_BAND.LT.IW)N_BAND=IW20 CONTINUEN_BAND=(N_BAND+1)*2IF(ID.EQ.1) THENELSEPE=PE/(1.0-PR*PR)PR=PR/(1.0-PR)END IFR ETURNENDcC to form the stiffness matrix of elementSUBROUTINE FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE) DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6),DD(3,3), & AKE(6,6), SS(6,6)CALL CAL_DD(PE,PR,DD)CALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)DO 10 I=1,3DO 10 J=1,6SS(I,J)=0.0DO 10 K=1,310 SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)DO 20 I=1,6DO 20 J=1,6AKE(I,J)=0.0DO 20 K=1,320 AKE(I,J)=AKE(I,J)+SS(K,I)*BB(K,J)*AE*PTRETURNENDcc to form banded global stiffness matrixSUBROUTINE BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,IJK_ELE,X,Y,PE, & PR,PT,AK)DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),AKE(6,6),AK(500,100)N_DOF=2*N_NODEDO 40 I=1,N_DOFDO 40 J=1,N_BAND40 AK(I,J)=0DO 50 IE=1,N_ELECALL FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE)DO 50 I=1,3DO 50 II=1,2IH=2*(I-1)+IIIDH=2*(IJK_ELE(IE,I)-1)+IIDO 50 J=1,3DO 50 JJ=1,2IL=2*(J-1)+JJIZL=2*(IJK_ELE(IE,J)-1)+JJIDL=IZL-IDH+1IF(IDL.LE.0) THENELSEAK(IDH,IDL)=AK(IDH,IDL)+AKE(IH,IL)END IF50 CONTINUERETURNENDcc to calculate the area of elementSUBROUTINE CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE)I=IJK_ELE(IE,1)J=IJK_ELE(IE,2)K=IJK_ELE(IE,3)XIJ=X(J)-X(I)YIJ=Y(J)-Y(I)XIK=X(K)-X(I)YIK=Y(K)-Y(I)AE=(XIJ*YIK-XIK*YIJ)/2.0RETURNENDcc to calculate the elastic matrix of elementSUBROUTINE CAL_DD(PE,PR,DD)DIMENSION DD(3,3)DO 10 I=1,3DO 10 J=1,310 DD(I,J)=0.0DD(1,1)=PE/(1.0-PR*PR)DD(1,2)=PE*PR/(1.0-PR*PR)DD(2,1)=DD(1,2)DD(2,2)=DD(1,1)DD(3,3)=PE/((1.0+PR)*2.0)RETURNENDcc to calculate the strain-displacement matrix of elementSUBROUTINE CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB) DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6)I=IJK_ELE(IE,1)J=IJK_ELE(IE,2)K=IJK_ELE(IE,3)DO 10 II=1,3DO 10 JJ=1,310 BB(II,JJ)=0.0BB(1,1)=Y(J)-Y(K)BB(1,3)=Y(K)-Y(I)BB(1,5)=Y(I)-Y(J)BB(2,2)=X(K)-X(J)BB(2,4)=X(I)-X(K)BB(2,6)=X(J)-X(I)BB(3,1)=BB(2,2)BB(3,2)=BB(1,1)BB(3,3)=BB(2,4)BB(3,4)=BB(1,3)BB(3,5)=BB(2,6)BB(3,6)=BB(1,5)CALL CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)DO 20 I1=1,3DO 20 J1=1,620 BB(I1,J1)=BB(I1,J1)/(2.0*AE)RETURNENDcc to form the global load matrixSUBROUTINE FORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK, & RESULT_N)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),P_IJK(N_LOAD,3), & RESULT_N(N_DOF)DO 10 I=1,N_DOF10 RESULT_N(I)=0.0DO 20 I=1,N_LOADII=P_IJK(I,1)RESULT_N(2*II-1)=P_IJK(I,2)20 RESULT_N(2*II)=P_IJK(I,3)RETURNENDcc to deal with BC(u) (here only for fixed displacement) using "1-0" method SUBROUTINE DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N) DIMENSION RESULT_N(N_DOF),IJK_U(N_BC,3),AK(500,100)DO 30 I=1,N_BCIR=IJK_U(I,1)DO 30 J=2,3IF(IJK_U(I,J).EQ.0)THENELSEII=2*IR+J-3AK(II,1)=1.0RESULT_N(II)=0.0DO 10 JJ=2,N_BAND10 AK(II,JJ)=0.0DO 20 JJ=2,II20 AK(II-JJ+1,JJ)=0.0END IF30 CONTINUERETURNENDcc to solve the banded FEM equation by GAUSS eliminationSUBROUTINE SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N) DIMENSION RESULT_N(N_DOF),AK(500,100)DO 20 K=1,N_DOF-1IF(N_DOF.GT.K+N_BAND-1)IM=K+N_BAND-1IF(N_DOF.LE.K+N_BAND-1)IM=N_DOFDO 20 I=K+1,IML=I-K+1C=AK(K,L)/AK(K,1)IW=N_BAND-L+1DO 10 J=1,IWM=J+I-K10 AK(I,J)=AK(I,J)-C*AK(K,M)20 RESULT_N(I)=RESULT_N(I)-C*RESULT_N(K)RESULT_N(N_DOF)=RESULT_N(N_DOF)/AK(N_DOF,1)DO 40 I1=1,N_DOF-1I=N_DOF-I1IF(N_BAND.GT.N_DOF-I-1)JQ=N_DOF-I+1IF(N_BAND.LE.N_DOF-I-1)JQ=N_BANDDO 30 J=2,JQK=J+I-130 RESULT_N(I)=RESULT_N(I)-AK(I,J)*RESULT_N(K)40 RESULT_N(I)=RESULT_N(I)/AK(I,1)WRITE(8,50)50 FORMAT(/12X,'* * * * * RESULTS BY FEM2D * * * * *',//8X,&'--DISPLACEMENT OF NODE--'//5X,'NODE NO',8X,'X-DISP',8X,'Y-DISP') DO 60 I=1,N_NODE60 WRITE(8,70) I,RESULT_N(2*I-1),RESULT_N(2*I)70 FORMAT(8X,I5,7X,2E15.6)RETURNENDcc calculate the stress components of element and nodeSUBROUTINECAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N, &STS_ELE,STS_ND)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),DD(3,3),BB(3,6), &SS(3,6),RESULT_N(N_DOF),DISP_E(6)DIMENSION STS_ELE(500,3),STS_ND(500,3)WRITE(8,10)10 FORMAT(//8X,'--STRESSES OF ELEMENT--')CALL CAL_DD(PE,PR,DD)DO 50 IE=1,N_ELECALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)DO 20 I=1,3DO 20 J=1,6SS(I,J)=0.0DO 20 K=1,320 SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)DO 30 I=1,3DO 30 J=1,2IH=2*(I-1)+JIW=2*(IJK_ELE(IE,I)-1)+J30 DISP_E(IH)=RESULT_N(IW)STX=0STY=0TXY=0DO 40 J=1,6STX=STX+SS(1,J)*DISP_E(J)STY=STY+SS(2,J)*DISP_E(J)40 TXY=TXY+SS(3,J)*DISP_E(J)STS_ELE(IE,1)=STXSTS_ELE(IE,2)=STYSTS_ELE(IE,3)=TXY50 WRITE(8,60)IE,STX,STY,TXY60 FORMAT(1X,'ELEMENT NO.=',I5/18X,'STX=',E12.6,5X,'STY=',&E12.6,2X,'TXY=',E12.6)c the following part is to calculate stress components of nodeWRITE(8,55)55 FORMAT(//8X,'--STRESSES OF NODE--')DO 90 I=1,N_NODEA=0.B=0.C=0.II=0DO 70 K=1,N_ELEDO 70 J=1,3IF(IJK_ELE(K,J).EQ.I) THENII=II+1A=A+STS_ELE(K,1)B=B+STS_ELE(K,2)C=C+STS_ELE(K,3)END IF70 CONTINUESTS_ND(I,1)=A/IISTS_ND(I,2)=B/IISTS_ND(I,3)=C/IIWRITE(8,75)I,STS_ND(I,1),STS_ND(I,2),STS_ND(I,3)75 FORMAT(1X,'NODE NO.=',I5/18X,'STX=',E12.6,5X,'STY=',&E12.6,2X,'TXY=',E12.6)90 CONTINUERETURNENDc FEM2D programm end[算例结果]：chengxu.for所输出的数据文件DATA.OUT数据内容如下：=========PLANE STRESS PROBLEM========* * * * * RESULTS BY FEM2D * * * * *--DISPLACEMENT OF NODE--NODE NO X-DISP Y-DISP1 .000000E+00 -.525275E+012 .000000E+00 -.225275E+013 -.108791E+01 -.137363E+014 .000000E+00 .000000E+005 -.824176E+00 .000000E+006 -.182418E+01 .000000E+00--STRESSES OF ELEMENT--ELEMENT NO.= 1STX=-.108791E+01 STY=-.300000E+01 TXY= .439560E+00ELEMENT NO.= 2STX=-.824176E+00 STY=-.225275E+01 TXY= .000000E+00ELEMENT NO.= 3STX=-.108791E+01 STY=-.137363E+01 TXY= .307692E+00ELEMENT NO.= 4STX=-.100000E+01 STY=-.137363E+01 TXY=-.131868E+00--STRESSES OF NODE--NODE NO.= 1STX=-.108791E+01 STY=-.300000E+01 TXY= .439560E+00NODE NO.= 2STX=-.100000E+01 STY=-.220879E+01 TXY= .249084E+00NODE NO.= 3STX=-.105861E+01 STY=-.191575E+01 TXY= .205128E+00NODE NO.= 4STX=-.824176E+00 STY=-.225275E+01 TXY= .000000E+00NODE NO.= 5STX=-.970696E+00 STY=-.166667E+01 TXY= .586081E-01NODE NO.= 6STX=-.100000E+01 STY=-.137363E+01 TXY=-.131868E+00[结论与体会]:通过本次的课程设计，我对有限元的概念有了更加深刻的理解，同时也弥补了平时学习是疏忽的地方，充实了有限元知识。

MATLAB期末⼤作业学号：姓名：《Matlab/Simulink在数学计算与仿真中的应⽤》⼤作业1.假设地球和⽕星绕太阳运转的半径分别为r和2r，利⽤comet指令动画显⽰从地球到⽕星的转移轨迹（r可以任意取值，要求实时显⽰探测器、太阳、地球和⽕星的位置）。

解函数function comet(varargin)[ax,args,nargs] = axescheck(varargin{:});error(nargchk(1,3,nargs,'struct'));% Parse the rest of the inputsif nargs < 2, x = args{1}; y = x; x = 1:length(y); endif nargs == 2, [x,y] = deal(args{:}); endif nargs < 3, p = 0.10; endif nargs == 3, [x,y,p] = deal(args{:}); endif ~isscalar(p) || ~isreal(p) || p < 0 || p >= 1error('MATLAB:comet:InvalidP', ...'The input ''p'' must be a real scalar between 0 and 1.'); End指令 %particle_motiont = 0:5:16013;r1=6.7e6;%随便给定参数%---------------------------r2=2*r1;g=9.8;R=6.378e6;m=g*R^2;%内轨道v_inner=sqrt(m/r1);w_inner=v_inner/r1;x_inter=r1*cos(w_inner*t);y_inter=r1*sin(w_inner*t);%外轨道v_outer=sqrt(m/r2);w_outer=v_outer/r2;x_outer=r2*cos(w_outer*t);y_outer=r2*sin(w_outer*t);%控制器转移轨道a=(r1+r2)/2;E=-m/(2*a);V_near=sqrt(m*(2/r1-2/(r1+r2)));%转移轨道在近地点的速度V_far=sqrt(m*(2/r2-2/(r1+r2)));%转移轨道在远地点的速度h=r1*V_near;%由于在近地点的速度垂直于位置失量, h是转移轨道的⽐动量矩e=sqrt(1+2*E*h^2/m^2);%e为椭圆轨迹的偏⼼率TOF=pi*sqrt(a^3/m);%转移轨道是椭圆的⼀半及飞⾏时间是周期的⼀半(开普勒第三定律) w=pi/TOF;%椭圆轨迹的⾓速度c=a*e;b=sqrt(a^2-c^2);x_ellipse=a*cos(w*t)-0.5*r1;y_ellipse=b*sin(w*t);%动画显⽰运动轨迹x=[x_inter;x_outer;x_ellipse]';y=[y_inter;y_outer;y_ellipse]';comet(x,y)%---------------------------动态图像如下：2．利⽤两种不同途径求解边值问题dfdxf gdgdxf g f g=+=-+==34430001,,(),()的解.途径⼀：指令syms f g[f,g]=dsolve('Df=3*f+4*g,Dg=-4*f+3*g','f(0)=0,g(0)=1'); disp('f=');disp(f)disp('g=');disp(g)结果（Matlab 7.8版本）f=i/(2*exp(t*(4*i - 3))) - (i*exp(t*(4*i + 3)))/2g=exp(t*(4*i + 3))/2 + 1/(2*exp(t*(4*i - 3)))（Matlab 6.5版本）f=exp(3*t)*sin(4*t)g=exp(3*t)*cos(4*t)>>途径⼆：%problem2function dy=problem2(t,y)dy = zeros(2,1);dy(1) = 3*y(1)+4*y(2);dy(2) = -4*y(1)+3*y(2);[t,y] = ode45('problem2',[0 2],[0 1]);plot(t,y(:,1),'r',t,y(:,2),'b');图23．假设著名的Lorenz 模型的状态⽅程表⽰为-+-=+-=+-=)()()()()()()()()()()()(322133223211t x t x t x t x t x t x t x t x t x t x t x t x σρρβ其中，设28,10,3/8===σρβ。

MATLAB作业⼀、必答题：1. MATLAB系统由那些部分组成？答：MATLAB系统主要由开发环境、MATLAB语⾔、MATLAB数学函数库、图形功能和应⽤程序接⼝五个部分组成。

2. 如何启动M⽂件编辑/调试器？答：在操作界⾯上选择“建⽴新⽂件”或“打开⽂件”操作时，M⽂件编辑/调试器将被启动。

在命令窗⼝中键⼊“edit”命令也可以启动M⽂件编辑/调试器。

3. 存储在⼯作空间中的数组能编辑吗？如何操作？答：存储在⼯作空间的数组可以通过数组编辑器进⾏编辑：在⼯作空间浏览器中双击要编辑的数组名打开数组编辑器，再选中要修改的数据单元，输⼊修改内容即可。

4. 在MATLAB中有⼏种获得帮助的途径？答：在MATLAB中有多种获得帮助的途径：（1）帮助浏览器：选择view菜单中的Help菜单项或选择Help菜单中的MATLAB Help菜单项可以打开帮助浏览器；（2）help命令：在命令窗⼝键⼊“help” 命令可以列出帮助主题，键⼊“help 函数名”可以得到指定函数的在线帮助信息；（3）lookfor命令：在命令窗⼝键⼊“lookfor 关键词”可以搜索出⼀系列与给定关键词相关的命令和函数（4）模糊查询：输⼊命令的前⼏个字母，然后按Tab键，就可以列出所有以这⼏个字母开始的命令和函数。

5. 有⼏种建⽴矩阵的⽅法？各有什么优点？答：（1）以直接列出元素的形式输⼊；（2）通过语句和函数产⽣；（3）.在m⽂件中创建矩阵；（4）从外部的数据⽂件中装⼊。

6. 命令⽂件与函数⽂件的主要区别是什么？答：命令⽂件: M⽂件中最简单的⼀种，不需输出输⼊参数，⽤M ⽂件可以控制⼯作空间的所有数据。

运⾏过程中产⽣的变量都是全局变量。

运⾏⼀个命令⽂件等价于从命令窗⼝中顺序运⾏⽂件⾥的命令，程序不需要预先定义，只要依次将命令编辑在命令⽂件中即可。

函数⽂件：如果M⽂件的第⼀个可执⾏⾏以function开始，便是函数⽂件，每⼀个函数⽂件定义⼀个函数。