matlab 微积分基本运算
§1 解方程和方程组解
1. 线性方程组求解
对于方程 AX = B ,其中 A 是( m ×n )的矩阵有三种情形： 1)当n=m 且A 非奇异时,此方程为“恰定”方程组。 2)当 n > m 时,此方程为“超定”方程组。 3)当n 2)采用左除运算解方程 x=A\B 例1 “求逆”法和“左除”法求下列方程组的解 ???? ? ????=+=++=++=++=+1 5065065065165545434323212 1x x x x x x x x x x x x x 在Matlab 编辑器中建立M 文件fanex1.m ： A=[5 6 0 0 0 1 5 6 0 0 0 1 5 6 0 0 0 1 5 6 0 0 0 1 5]; B=[1 0 0 0 1]'; R_A=rank(A) %求秩 X1=A\B %用"左除"法解恰定方程所得的解 X2=inv(A)*B %用"求逆"法解恰定方程所得的解 运行后结果如下 R_A = 5 X1 = 2.2662 -1.7218 1.0571 -0.5940 0.3188 X2 = 2.2662 -1.7218 1.0571 -0.5940 0.3188 两种方法所求方程组的解相同。 (2)MATLAB 解超定方程AX=B 的方法 对于方程 AX = B ,其中 A 是( m ×n )的矩阵， n > m ，如果A 列满秩，则此方程是没有精确解的。然而在实际工程应用中，求得其最小二乘解也是有意义的。基本解法有： 1)采用求伪逆运算解方程 x=pinv(A)*B 说明：此解为最小二乘解x=inv(A ’*A)*A*B,这里pinv(A) =inv(A ’*A)*A. 2)采用左除运算解方程 x=A\B 例2 “求伪逆”法和“左除”法求下列方程组的解 ??? ??=+=+=+1 221421 221 2121x x x x x x 命令如下： >> a=[1 2;2 4;2 2]; >> b=[1,1,1]'; >> xc=a\b %用左除运算解方程 运行得结果: xc = 0.4000 0.1000 >> xd=pinv(a)*b %用求伪逆运算解方程 运行得结果: xd = 0.4000 0.1000 >> a*xc-b %xc 是否满足方程ax=b 运行得结果: ans = -0.4000 0.2000 0.0000 可见xc 并不是方程的精确解。 (3) MATLAB 解欠定方程AX=B 的方法 欠定方程从理论上说是有无穷多个解的，如果利用求“伪逆”法和“左除”法来求解，只能得到其中一个解。基本方法： 1)采用求伪逆运算解方程 x=pinv(A)*B 2)采用左除运算解方程 x=A\B 例3 求方程组 ??? ??=--+=+--=--+0 8954433134321 43214321x x x x x x x x x x x x 的一个特解 解：在Matlab 编辑器中建立M 文件：fanex2.m A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8]; B=[1 4 0]'; X=A\B X = 0 0 -0.5333 0.6000 (4) MATLAB 求线性齐次方程组的通解 在Matlab 中，函数null 用来求解零空间，即满足AX=0的解空间，实际上是求出解空间的一组基（基础解系）。基本格式： 1)格式：z = null(A) 说明： z 的列向量为方程组的正交规范基，满足I Z Z =?' 2)格式：z = null(A ，’r ’) 说明: z 的列向量是方程AX=0的有理基. 例4 求解方程组的通解： ??? ??=---=--+=+++0 3402220224321 43214321x x x x x x x x x x x x 解：在Matlab 编辑器中建立M 文件：fanexm A=[1 2 2 1;2 1 -2 -2;1 -1 -4 -3]; format rat %指定有理式格式输出 B=null(A,'r') %求解空间的有理基 运行后显示结果如下： B = 2 5/ 3 -2 -4/3 1 0 0 1 写出通解： syms k1 k2 X=k1*B(:,1)+k2*B(:,2) %写出方程组的通解 pretty(X) %让通解表达式更加精美 运行后结果如下： X = [ 2*k1+5/3*k2] [ -2*k1-4/3*k2] [ k1] [ k2] % 下面是其简化形式 [2k1 + 5/3k2 ] [ ] [-2k1 - 4/3k2] [ ] [ k1 ] [ ] [ k2 ] (4)MATLAB 求非齐次线性方程组的通解 非齐次线性方程组需要先判断方程组是否有解，若有解，再去求通解。因此，步骤为： 第一步：判断AX=b 是否有解，若有解则进行第二步 第二步：求AX=b 的一个特解 第三步：求AX=0的通解 第四步：AX=b 的通解= AX=0的通解+AX=b 的一个特解。 例5 求解方程组 ??? ??=-++=-+-=-+-3 22223531324321 43214321x x x x x x x x x x x x 解：在Matlab 编辑器中建立M 文件：fanex5.m A=[1 -2 3 -1;3 -1 5 -3;2 1 2 -2]; b=[1 2 3]'; B=[A b]; n=4; R_A=rank(A) R_B=rank(B) format rat if R_A==R_B&R_A==n %判断有唯一解 X=A\b elseif R_A==R_B&R_A C=null(A,'r') %求AX=0的基础解系 else X='equition no solve' %判断无解 end 运行后结果显示： R_A = 2 R_B = 3 X = equition no solve 说明该方程组无解 例6 求解方程组的通解： ??? ??=--+=+--=--+0 8954433134321 43214321x x x x x x x x x x x x 解法一：在Matlab 编辑器中建立M 文件：fanex6.m A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8]; b=[1 4 0]'; B=[A b]; n=4; R_A=rank(A) R_B=rank(B) format rat if R_A==R_B&R_A==n X=A\b elseif R_A==R_B&R_A C=null(A,'r') else X='Equation has no solves' end 运行后结果显示为： R_A = 2 R_B = 2 X = 0 0 -8/15 3/5 C = 3/2 -3/4 3/2 7/4 1 0 0 1 所以原方程组的通解为 X=k 1??????? ??012/32/3+k 2 ??????? ??-104/74/3+??????? ? ?-5/315/800 解法二：在Matlab 编辑器中建立M 文件：LX07212.m A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8]; b=[1 4 0]'; B=[A b]; C=rref(B) %求增广矩阵的行最简形，可得最简同解方程组。 运行后结果显示为： C = 1 0 -3/ 2 3/4 5/4 0 1 -3/2 -7/4 -1/4 0 0 0 0 0 对应齐次方程组的基础解系为： ??? ???? ??=012/32/31ξ， ?????? ? ??-=104/74/32ξ 非齐次方程组的特解为： ?????? ? ??-=004/14/5*η 所以，原方程组的通解为： X=k 11ξ+k 22ξ+*η 2.代数方程组求解 在MMlab 中解方程组的命令有linsolve 和solve ，其中命令linsolve 专门用于求解线性方程组，而命令solve 可适用于所有的代数方程(组)。下面分别加以介绍。 (1) 求解线性方程组linsolve linsolve 命令对解形如AX=B 的线性方程组运用自如。但在对矩阵A 的运算中有如下限制：A 必须至少是行满秩的。Linsolve 具体使用格式为： X=linsolve(A,B) 功能：得到线性方程组AX=B 的特解X 。 例7 求线性方程组??? ??=+-=++=+-20 141861234834483321 321321x x x x x x x x x AX=b 的解。 命令如下： >>A=[3,-8,4;3,48,3;6,-18,14]; >>b=[4;12;20]; >>X=linsolve(A,b) %调用linsolve 函数求特解 运行后结果显示为： X = [ -156/167] [ 30/167] [ 344/167] >>A\b %用另一种方法求特解 运行后结果显示为： ans = -0.9341 0.1796 2.0599 例8 求解方程组的特解： ??? ??=--+=+--=--+0 8954433134321 43214321x x x x x x x x x x x x 命令如下： >>A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8]; >>b=[1 4 0]'; >>X=linsolve(A,b) 运行后结果显示为： Warning: System is rank deficient. Solution is not unique. > In E:\MATLAB6p5\toolbox\symbolic\@sym\mldivide.m at line 38 In E:\MATLAB6p5\toolbox\symbolic\linsolve.m at line 8 X = [ 5/4] [ -1/4] [ 0] [ 0] (2) 一般符号表达式表示的代数方程求解 在MATLAB 中，求解用符号表达式表示的代数方程可由函数solve 实现，其调用格式为： 1) 格式：x=solve(S) 功能：求解符号表达式S 的代数方程，求解变量为默认变量v=findsym(S),将结果赋给x 。其中S 为包含方程(一个)等式的字符串(可以是函数名，或者是描述方程的字符串，关于方程组的输入在前面章节中有讨论)； 2) 格式：x=solve(S,v) 功能：求解符号表达式S 的代数方程，求解变量为v, 将结果赋给x 。 3) 格式：[x1,x2, …,xn]=solve(S1,S 2,…,Sn) 功能：求解符号表达式S1, S 2,…, Sn 组成的代数方程组，对n 个默认变量求解，将结果赋给[x1,x2, …,xn]。 4) 格式：solve(S1,S 2,…,S n,v1,v2,…,vn) 功能：求解符号表达式S1, S 2,…, Sn 组成的代数方程组，求解变量分别v1,v2,…,vn，将结果赋给[x1,x2, …,xn]。 例9 解方程: (1) 2 2 144212-+ =-++x x x x (2)1743 3=--- x x x ; (3)1)4 3sin(2=- π x ; (4)10=+x xe x ;。 命令如下： >> X1=solve('1/(x+2)+4*x/(x^2-4)=1+2/(x-2)','x') %解方程(1) 运行后结果显示为： X1 = 1 >> f=sym('x-(x^3-4*x-7)^(1/3)=1'); >> x2=solve(f) %解方程(2) 运行后结果显示为： x2 = 3 >> x3=solve('2*sin(3*x-pi/4)=1') %解方程(3) 运行后结果显示为： x3 = 5/36*pi >> x4=solve('x+x*exp(x)-10','x') %解方程(4)。 运行后结果显示为： x4 = 1.6335 例10 解方程组 ???=-+=-+0 60 622x y y x 命令如下： >> clear >> [x,y]=solve('x^2+y-6=0','y^2+x-6=0','x','y') 运行后结果显示为： x = [ 2] [ -3] [ 1/2-1/2*21^(1/2)] [ 1/2+1/2*21^(1/2)] y = [ 2] [ -3] [ 1/2+1/2*21^(1/2)] [ 1/2-1/2*21^(1/2)] §2 数列极限、函数极限及连续 本部分内容主要讲述极限、无穷大、无穷小、连续等概念和极限的四则运算、 数列和函数极限的常用求法。 1. 求极限 极限问题的类型大概有数列的极限（n n a ∞ →lim ）、一元函数的极限(包括无穷极限：)(lim ),(lim ),(lim x f x f x f x x x +∞ →-∞ →∞ →和 有限点极限：)(lim 0 x f x x →及左、右极限)(lim ),(lim 0 x f x f x x x x +→-→)和二元(多元)函数的极限(),(lim 0 y x f y y x x →→)。以上是数学 中出现的不同形式的极限。但求极限的命令在Matlab 中是相同的，都是由命令limit 来完成。该函数的格式及功能有： (1) 格式：limit(f, 'v',a) 功能：求符号函数f 的极限值。即计算当变量v 趋近于常数a 时，f 函数的极限值。 (2) 格式：limit(f) 功能：求符号函数f 的极限值。符号函数f 的变量为函数findsym(f)确定的默认变量；没有指定变量的目标值时，系统默认变量趋近于0，即a=0的情况。 (3) 格式： limit(f,a) 功能：求符号函数f 的极限值。由于没有指定符号函数f 的自变量，则使用该格式时，符号函数f 的变量为函数findsym(f)确定的默认自变量，既变量'v'趋近于a 。 (4) 格式：limit(f, 'v',a,'right') 功能：求符号函数f 的极限值。'right'表示变量'v 从右边趋近于a 。 (5) 格式：limit(f, 'v',a,'left') 功能：求符号函数f 的极限值。'left'表示变量v 从左边趋近于a 。 例1 计算下列极限 (1) 116lim 232 +++-∞→n n n n n ； (2) a x a x m m a x --→11 lim ； (3) ))1(1 ln 1( lim 2 21 --→x x x x ； (4) 2 2 lim a x a x a x a x ----+ →； (5)x x x a 3)21(lim + ∞ →； (6) x x x sin lim 0→； (7) 2 2 20 1 3) cos ln(lim y x y x y x ++→→ (2 2 21 02 2 20 13)cos ln(lim (lim 3)cos ln(lim y x y x y x y x x y y x ++=++→→→→) Matlab 程序如下： >> syms x y n a m; %定义x 、n 、 a 、 m 为符号变量，其中分号的作用是不显示结果. >> y1=(6*n^2-n+1)/(n^3-n^2+2); %定义符号表达式 >> Lim_y1=limit(f,n,inf) %求1 1 6lim 232+++-∞→n n n n n 极限 结果为： Lim_y1= 0 >> Lim_y2= limit( (x^(1/m)-a^(1/m))/(x-a),x,a) %求 a x a x m m --11极限 结果为： Lim_y2 = a^(1/m)/a/m >> Lim_y3=limit(1/(x*log(x)^2)-1/(x-1)^2,1) %求)) 1(1 ln 1( lim 2 21 --→x x x x 极限 结果为： Lim_y3 = 1/12 >>y4=(sqrt(x)-sqrt(a)-sqrt(x-a))/sqrt(x*x-a*a); %定义符号表达式 >> Lim_y4=limit(y4,x,a,'right') %求2 2 lim a x a x a x a x ----→极限 结果为： Lim_y4= -1/2*2^(1/2)/a^(1/2) >> Lim_y5=limit ((1+2*a/x)^(3*x) , x, inf ) %求x x x a 3)21(lim + ∞→极限 结果为： Lim_y5= exp(6*a) >> Lim_y6=limit(sin(x)/x) %求x x x sin lim 0→极限 结果为： Lim_y6 = 1 >> Lim_y7=limit(limit((log(x^2+cos(y))/sqrt(3*x^2+y^2),x,1),y,0) 结果为： Lim_y7 = 1/3*log(2)*3^(1/2) 注意：求极限的基本步骤为： 1、把数学形式的表达式转化成Matlab 的表达式。当然表达式由多个符号组成．并且符号都不是代表某一具体的数值，其值是一个变动的量，即所谓的符号变量。为此，我们要建立一些必要的符号对象，如：符号变量、符号表达式等。 2、使用求极限的命令limit 来对符号表达式进行计算，其中要指明是对哪一个符号变量做极限运算；若不指定，则符号变量由命令findsym(f)指定。 3、要指定符号变量的“极限点”以及趋于“极限点”的方向；若不指定“极限点”的值，则缺省地认为是0点。 4、多元函数的极限要比一元函数的极限复杂，因为按照不同的方式趋近于某一目标，其结果可能不一样，如果结果不同则极限就不存在。基于这样的事实，没有提供专门的命令函数来计算多元函数的极限，这一功能仍然是由命令函数limit 来完成，但只能求在理论上已经证明积分与路径无关的函数，即把二重极限化成二次极限来计算。 2. 连续性讨论 函数)(x f 在点a x =连续，定义为： )()(lim a f x f a x =→(或)()(lim )(lim a f x f x f a x a x ==-+→→) 即函数)(x f 在点a x =处极限值等于函数值(或左、右极限存在并相等，而且等于函数值)。 例2 讨论下列函数在指定点的连续性： (1) 函数?????-=--≠+--=171,1 65)(2x x x x x x f 在1-=x 处的连续性； (2) 函数?? ? ??<≥--=0sin 0,51)(2x x x x x x x g 在0=x 处的连续性； 在MATLAB 命令窗口，输入命令： >> syms x; %定义x 为符号变量 >> f=(x^2-5*x-6)/(x+1); %定义符号表达式 >>Lim_f=limit (f, x, -1) %求极限 结果为： Lim_f = -7 显然极限值等于函数值，故)(x f 在在1-=x 处的连续。 >> left_lim_g=limit(sin(x)/x , x, 0, ’left’) 结果为： left_lim_g = 1 >> right_lim_g=limit(1-x^2-5*x ,x, 0, ’right’) 结果为： right_lim_g = 1 左、右极限存在并相等，而且等于函数值, 故连续。 3．简单应用 例3细菌繁殖问题 问题 由实验知，某种细菌繁殖的速度在培养基充足等条件满足时与当时已有的数量A0成正比，即V=kA0（k ＞0为比例常数），问经过时间t 以后细菌的数量是多少？ 问题求解： 为了计算出t 时的数量，我们将时间间隔【0，t 】分为n 等分。由于细菌的繁殖是连续变化的，在很短的一段时间内数量的变化很小，繁殖速度可近似看作不变。因此，在第一段时间【0，t/n 】末细菌繁殖的数量为：kA0t/n ，第 一段时间末细菌的数量为：A0（1+kt/n ）；同样，第二段时间末细菌的数量为：A0(1+k n t )2，…依次类推，到最后一段时间末细菌的数量为A0(1+k n t )n 。 显然，这是一个近似值，因为我们假设了在每一小段时间【（i-1）t/n ，it/n 】（i ＝1，2，3…，n ）内细菌繁殖的速度不变（同时还假设了各小段时间内只繁殖一次）。可以看出，当时间间隔分得越细（即n 越大）时这个值越接近精确值，若对时间间隔无限细分（即n →∞），则可求得其精确值。所以经过时间t 后细菌的总数是 lim ∞ →n A0(1+k n t )n ＝A0lim ∞→n (1+k n t )n x 1n kt ＝令 A0lim ∞→x ［(1+x 1 )x ］kt ＝A0［ lim ∞ →x (1+x 1 )x ］kt ＝A0 ekt 事实上，现实世界中不少事物的生长规律都服从这个模型，所以也称y ＝Aekt 为生长函数。下面看一个实例。 例 已知一种细菌的个数按指数方式增长，收集到的数据是5天后的细菌个数为936，10天后细菌个数为2190。问： （1）开始时细菌个数是多少？ （2）如果继续以现在的速度增长下去，60天后细菌的个数是多少？ 解 细菌繁殖服从生长函数y ＝A0ekt 。由题目所给数据，得 ?????==10k 05k 0e A 2190e A 936 解此方程组，得A0＝400，k ＝0.17。即开始时细菌个数为400。按此速度增长下去，则60天后细菌个数为 Y （60）＝400*e60*0.17＝10761200. 注 这里仅两组数据去定A0，k 必有误差。为了得到较为准确的A0，k 的估计值A0,k ，应当多收集一些数据，然后用以后将介绍的最小二乘法确定之。 现由matlab 符号和数值求解如下 ： 由符号计算求解： syms A k t n; %定义符号变量 lim=limit(A*(1+k*t/n)^n,n,inf); %求极限 由此可得出经过时间t 后细菌总数为： lim=exp(k*t)*A 具体计算求解： 例 已知一种细菌的个数按指数方式增长，收集到的数据是5天后细菌个数为936， 10天后细菌个数为2190。问： (1)开始时细菌个数是多少? (2)如果继续以现在的速度增长下去,60天后细菌的个数是多少? 解 细菌繁殖服从生产函数y=Aekt.由题目所给数据得 936=Ae5k 2190=Ae10k 解此方程组即可求解,现用符号求解如下: syms A k; %定义符号变量A 和k [A k]=solve('936=A*exp(5*k)','2190=A*exp(10*k)','A,k'); %用符号方程求解A 和k A=double(A) %将符号表示方法转变为数值表示方法 k=double(k) y=A*exp(60*k); %求出60天后的细菌个数 y=double(y) %将符号表示方法转变为数值表示方法 例4 某储户将10万元的人民币以活期的形式存入银行，年利率为5%0，如果银行允许储户在一年内可任意次结算，在不计利息税的情况下，若储户等间隔地结算n 次，每次结算后将本息全部存入银行，问一年后该储户的本息和是多少？随着结算次数的无限增加，一年后该储户是否会成为百万富翁？ 问题分析： 若该储户每季度结算一次，则每季度利率为： 4005 .0 故第一季度后储户本息共计： )4005 .01(100000+ 第二季度后储户本息共计： 2)4005.01(100000+ … … ，一年后储户本息共计： 4 )4005.01(100000+ 若该储户每月结算一次，则每月利率为：12005 .0 按上面的方法知一年后储户本息共计： 12 )12005.01(100000+ 若该储户每5天结算一次，则一年后本息共计 ：73 )73005.01(100000+ 若该储户的间隔地结算n 次，则一年后本息共计 ： n n )005.01(100000+ 随着结算次数的无限增加，即在上式中,故一年后本息共计： n n n )005.01(100000lim + ∞ → >> syms n; >> Limit[100000*(1 + 0.005/n)^n, n-> Infinity] 结果为： ans = 100000*exp(1/200) §3 函数的导数与最值 本部分内容主要讲述(偏)导数及其应用和最值的常用求法。 1． 函数求导 函数的导数问题的类型大致有一元函数的导数（一阶导数)(x f '和高阶导数)() (x f n ）多元函数的导数(包括无穷 极限：一阶偏导数y f x f ????,等、高阶偏导数n n n n y f x f ????,等、混合偏导数 y x f ???2等、复合函数求导、参数函数求导及隐函数求导。以上是数学中出现的不同形式的导数。但求导数的命令在Matlab 中是相同的，都是由命令diff 来完成。该 函数的格式及功能有： (1) 格式：diff(f, 'v',n) 功能：diff 函数求函数f 对变量'v'的n 阶(偏)导数。 (2) 格式：diff(f, 'v') 功能：对表达式f 中指定符号变量v 计算f 的1阶(偏)导数。 (3) 格式：diff(f) 功能：对表达式f 中的符号变量v 计算f 的1阶(偏)导数，其中v=findsym(f)。 (4) 格式：diff(f, n) 功能：对表达式f 中的符号变量v 计算f 的n 阶(偏)导数，其中v=findsym(f)。 (5) 格式：梯度jacobian(f,v) 功能：求向量函数)),,,(),,,,(),,,,((21211211n m n n x x x f x x x f x x x f f =的各分量的梯度： ? ?? ?????????????????????????????????n m m m n n x f x f x f x f x f x f x f x f x f 2 1 2221212111 2． 函数求导编程实例 (1)求一元函数的导数. １）)(x f y =的一阶导数. 例1 求dx x d 2 sin 导数； 打开matlab 指令窗，输入指令: >>x = sym('x'); %定义x 为符号变量 >>diff(sin(x^2)) %求dx x d 2 sin 导数 得结果: ans = 2*cos(x^2)*x 利用matlab 命令diff 一次可以求出若干个函数的导数. 例2 求下列函数的导数: (1)x x y 2cos 2cos 12 +=； (2)x x y =2； (3)x y sin 43=；(4)x y sin ln ln 4=。 输入命令: >> syms x >> a=diff([cos(x^2)+2*cos(2*x), x^x,4^(sin(x)),log(log(sin(x)))]) 得结果: a = [ -2*sin(x^2)*x-4*sin(2*x), x^x*(log(x)+1), 4^sin(x)*cos(x)*log(4), cos(x)/sin(x)/log(sin(x))] >>dy1_dx=a(1) 得结果: dy1_dx = -2*sin(x^2)*x-4*sin(2*x) >>dy2_dx=a(2) 得结果: dy2_dx = x^x*(log(x)+1) >>dy3_dx=a(3) 得结果: dy3_dx = 4^sin(x)*cos(x)*log(4) >>dy4_dx=a(4) 得结果: dy4_dx = cos(x)/sin(x)/log(sin(x)) 由本例可以看出，matlab 函数是对矩阵或向量进行操作的，a(i)表示向量a 的第i 个分量. ２）)(x f y =的高阶导数 例3 设x x x f 4cos )(3=,求)() 100(x f . 输入指令: >> syms x >> diff(x^3*cos(x),x,100) 得结果: ans = -970200*sin(x)-29700*x*cos(x)+300*x^2*sin(x)+x^3*cos(x) 3）参数方程所确定的函数的导数. 设参数方程? ? ?==)()(t y y t x x 确定函数)(x f y =，则y 的导数)() (t x t y dx dy ''= . 例4 求下列由参数方程确定的函数的一阶导数和二阶导数： ? ? ?=-=t t y t t x cos ) sin 1( 输入指令: >>syms t >>x=t*(1-sin(t)); >>y=t*cos(t); >>df1=diff(y,t)/diff(x,t) 得结果: df1 = (cos(t)-t*sin(t))/(1-sin(t)-t*cos(t)) >>df2=diff(df1,t)/diff(x,t) 得结果: df2 = ((-2*sin(t)-t*cos(t))/(1-sin(t)-t*cos(t))-(cos(t)-t*sin(t))/(1-sin(t)-t*cos(t))^2*(-2*cos(t)+t*s in(t)))/(1-sin(t)-t*cos(t)) (2) 求多元函数的偏导数 １）一阶偏导数 例5 设z y y xy u cos 3 ++=,求u 的一阶偏导数. 方法一、输入命令: >>syms x y z >>du_dx=diff(x*y+y^3+y*cos(z), x) 得结果: du_dx = y >> du_dx=diff(x*y+y^3+y*cos(z), y) 得结果: du_dx = x+3*y^2+cos(z) >> du_dx=diff(x*y+y^3+y*cos(z),z) 得结果: du_dx = -y*sin(z) 方法二、输入命令: >>syms x y z >>jacobian(x*y+y^3+y*cos(z),[x,y,z]) 得结果: ans = [ y, x+3*y^2+cos(z), -y*sin(z)] 给出矩阵),,( z u y u x u ??????. 说明：使用jacobian 命令求偏导数更为方便. 例6 求???? ? ??+=????? ??=z y x ye z xy z y x f z y x f z y x f z y x f z sin cos ),,(),,(),,(),,(2321jacobian 矩阵。 输入命令: >>syms x y z >> jacobian([x*y*cos(z),y*exp(z),x*sqrt(y^2+sin(z))],[x,y,z]) 得结果: ans = [ y*cos(z), x*cos(z), -x*y*sin(z)] [ 0, exp(z), y*exp(z)] [(y^2+sin(z))^(1/2), x/(y^2+sin(z))^(1/2)*y, 1/2*x/(y^2+sin(z))^(1/2)*cos(z)] 即得jacobian 矩阵： ? ? ?? ? ? ? ???? ???+++-=???????? ????????????????????????????z y z x z y xy z y ye e z xy z x z y z f y f x f z f y f x f z f y f x f z z sin 2cos sin sin 0sin cos cos 222 33322211 1 2）求高阶偏导数. 例7 求函数的导数: (i))(222y e x x -?? ; (ii) 2,022)(==-???y x x e y y x ； (iii) ? ?? ? ??????y e x y x sin 2 22; (i) 求)(222y e x x -?? >>syms x y >>diff(x^2*exp(-y),x,2) 得结果: ans = 2*exp(-y) (ii) 求2,022 )(==-???y x x e y y x 方法一：先求出 )(22x e y y x -???，再代值2,022)(==-???y x x e y y x 。 >>syms x y >>diff(diff(y^2*exp(-x),x,1),y,1) 得结果: ans = -2*y*exp(-x) >>x=0,y=2; >>-2*y*exp(-x) 得结果: ans = -4 即得结果： 4)(2,022 -=???==-y x x e y y x 。 方法二：利用定义 2 )()(lim )(2 ,020222 ,022 -?? -??=??? ==-=-→==-y e y x e y x e y y x y x x x x y y x x 输入命令: >>diff(y^2*exp(-x),x,1) 得结果: ans = -y^2*exp(-x) 输入命令: >>limit((-y^2+4)/(y-2),y,2) 得结果: ans = -4 即得结果： 4)(2,022 -=???==-y x x e y y x 说明：从此题可以看出，Matlab 的符号微分，只能进行符号运算，不能进行数值计算，即不能直接代值。 (iii)求??? ? ??????y e x y x sin 2 2 2 输入命令: >>syms x y >> diff(diff(exp(x^2)*sin(y),x,2),y,1) 得结果: ans = 2*exp(x^2)*cos(y)+4*x^2*exp(x^2)*cos(y) (3)求隐函数的导数或偏导数 例8 设20)sin(=+-x ye xy 确定函数)(x f y =，求 dx dy 。 实验01讲评、参考答案 讲评 未交实验报告的同学名单 批改情况: 问题1: 不仔细,式子中出错。 问题2: 提交的过程不完整。 问题3: 使用语句尾分号(;)不当,提交的过程中不该显示的结果显示。 问题4: 截屏窗口没有调整大小。 附参考答案: 实验01 MATLAB 运算基础 (第2章 MATLAB 数据及其运算) 一、实验目的 1、 熟悉启动与退出MATLAB 的方法。 2、 熟悉MATLAB 命令窗口的组成。 3、 掌握建立矩阵的方法。 4、 掌握MATLAB 各种表达式的书写规则以及常用函数的使用。 二、实验内容 1、 数学表达式计算 先求下列表达式的值,然后显示MATLAB 工作空间的使用情况并保存全部变量。 1、1 计算三角函数 12 2sin 851z e =+(注意:度要转换成弧度,e 2如何给出) 示例:点击Command Window 窗口右上角的,将命令窗口提出来成悬浮窗口,适当调 整窗口大小。 命令窗口中的执行过程: 《MATLAB 软件》课内实验 王平 1、2 计算自然对数 221 ln(1)2z x x =++,其中2120.45 5i x +??=??-??(提示:clc 命令擦除命令窗口,clear 则清除工作空间中的所有变量,使用时注意区别,慎用clear 命令。 应用点乘方) 命令窗口中的执行过程: 1、3 求数学表达式的一组值 0.30.330.3sin(0.3)ln , 3.0, 2.9, ,2.9,3.022 a a e e a z a a --+=++=-- 提示:利用冒号表达式生成a 向量,求各点的函数值时用点乘运算。 命令窗口中的执行过程: 用递推公式计算定积分 实验目的： 1．充分理解不稳定的计算方法会造成误差的积累，在计算过程中会导致误差的迅速增加，从而使结果产生较大的误差。 2．在选择数值计算公式来进行近似计算时，应学会选用那些在计算过程中不会导致误差迅速增长的计算公式。 3．理解不稳定的计算公式造成误差积累的来源及具体过程； 4．掌握简单的matlab语言进行数值计算的方法。 实验题目： 对n=0,1,2,…,20，计算定积分： 实验原理： 由于y(n)= = – 在计算时有两种迭代方法，如下： 方法一： y(n)=– 5*y(n-1),n=1,2,3, (20) 取y(0)= = ln6-ln5 ≈ 0.182322 方法二： 利用递推公式：y(n-1)=-*y(n),n=20,19, (1) 而且，由 = * ≤≤* = 可取：y(20)≈*()≈0.008730. 实验容： 对算法一,程序代码如下： function [y,n]=funa() syms k n t; t=0.182322; n=0; y=zeros(1,20); y(1)=t; for k=2:20 y(k)=1/k-5*y(k-1); n=n+1; end y(1:6) y(7:11) 对算法二，程序代码如下： %计算定积分； %n--表示迭代次数； %y用来存储结果； function [y,n]=f(); syms k y_20; y=zeros(21,1); n=1; y_20=(1/105+1/126)/2; y(21)=y_20; for k=21:-1:2 y(k-1)=1/(5*(k-1))-y(k)/5; n=n+1; end 实验结果： 由于计算过程中，前11个数字太小，后9个数字比较大，造成前面几个数字只显示0.0000的现象，所以先输出前6个，再输出7—11个，这样就能全部显示出来了。 算法一结果： [y,n]=funa %先显示一y(1)—y(6) ans = 0.1823 -0.4116 2.3914 -11.7069 58.7346 电子一班王申江 实验九数值微积分与方程数值求解 一、实验目的 1、掌握求数值导数和数值积分的方法 2、掌握代数方程数值求解的方法 3、掌握常微分方程数值求解的方法 二、实验内容 1、求函数在指定点的数值导数。 () 23 2 123,1,2,3 026 x x x f x x x x x == >>syms x >>f=[x x^2 x^3;1 2*x 3*x^2;0 2 6*x]; >>F=det(f) F=2*x^3 >>h=0.1 >>x=[0:h:4]; >>f=2*x^3; >>[dy,dx]=diff_ctr(f,h,1); >>y1=dy(dx==1) y1=6.0000 >>y2=dy(dx==2) y2=24.0000 >>y3=dy(dx==3) y3=54.0000 2、用数值方法求定积分。 （1） 210I π =?的近似值 a=inline('sqrt(cos(t.^2)+4*sin((2*t).^2)+1)'); I=quadl(a,0,2*pi) I = 6.7992 + 3.1526i （2）()1 202ln 11x I dx x +=+? b=inline('log(1+x)./(1+x.^2)'); I=quadl(b,0,1) I = 0.2722 3、分别用3种不同的数值方法解线性方程组。 6525494133422139211 x y z u x y z u x y z u x y u +-+=-??-+-=??++-=??-+=? A=[6,5,-2,5;9,-1,4,-1;3,4,2,-2;3,-9,0,2]; b=[-4,13,1,11]'; x=A\b 实验一 基本信号在MATLAB 中的表示和运算 一、实验目的 1． 学会用MA TLAB 表示常用连续信号的方法； 2． 学会用MA TLAB 进行信号基本运算的方法； 二、实验原理 1． 连续信号的MATLAB 表示 MATLAB 提供了大量的生成基本信号的函数，例如指数信号、正余弦信号。表示连续 时间信号有两种方法，一是数值法，二是符号法。数值法是定义某一时间范围和取样时间间 隔，然后调用该函数计算这些点的函数值，得到两组数值矢量，可用绘图语句画出其波形； 符号法是利用MATLAB 的符号运算功能，需定义符号变量和符号函数，运算结果是符号表 达的解析式，也可用绘图语句画出其波形图。 例1-1指数信号 指数信号在MATLAB 中用exp 函数表示。 如at Ae t f =)(，调用格式为 ft=A*exp(a*t) 程序是 A=1; a=-0.4; t=0:0.01:10; %定义时间点 ft=A*exp(a*t); %计算这些点的函数值 plot(t,ft); %画图命令，用直线段连接函数值表示曲线 grid on; %在图上画方格 例1-2 正弦信号 正弦信号在MATLAB 中用 sin 函数表示。 调用格式为 ft=A*sin(w*t+phi) A=1; w=2*pi; phi=pi/6; t=0:0.01:8; %定义时间点 ft=A*sin(w*t+phi); %计算这些点的函数值 plot(t,ft); %画图命令 grid on; %在图上画方格 例1-3 抽样信号 抽样信号Sa(t)=sin(t)/t 在MA TLAB 中用 sinc 函数表示。定义为 )/(sin )(πt c t Sa = t=-3*pi:pi/100:3*pi; ft=sinc(t/pi); plot(t,ft); grid on; axis([-10,10,-0.5,1.2]); %定义画图范围，横轴，纵轴 title('抽样信号') %定义图的标题名字 例1-4 三角信号 三角信号在MATLAB 中用 tripuls 函数表示。 一、符号积分 符号积分由函数int来实现。该函数的一般调用格式为： int(s)：没有指定积分变量和积分阶数时，系统按findsym函数指示的默认变量对被积函数或符号表达式s求不定积分； int(s,v)：以v为自变量，对被积函数或符号表达式s求不定积分； int(s,v,a,b)：求定积分运算。a,b分别表示定积分的下限和上限。该函数求被积函数在区间[a,b]上的定积分。a和b可以是两个具体的数，也可以是一个符号表达式，还可以是无穷(inf)。当函数f关于变量x在闭区间[a,b]上可积时，函数返回一个定积分结果。当a,b中有一个是inf时，函数返回一个广义积分。当a,b中有一个符号表达式时，函数返回一个符号函数。 例： 求函数x^2+y^2+z^2的三重积分。内积分上下限都是函数，对z积分下限是sqrt(x*y)，积分上限是x^2*y；对y积分下限是sqrt(x)，积分上限是x^2；对x的积分下限1，上限是2，求解如下： >>syms x y z %定义符号变量 >>F2=int(int(int(x^2+y^2+z^2,z,sqrt(x*y),x^2*y),y,sqrt(x),x^2),x,1,2) %注意定积分的书写格式 F2 = 1610027357/6563700-6072064/348075*2^(1/2)+14912/4641*2^(1/4)+64/225*2 ^(3/4) %给出有理数解 >>VF2=vpa(F2) %给出默认精度的数值解 VF2 = 224.92153573331143159790710032805 二、数值积分 1.数值积分基本原理 求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson)?法、牛顿－柯特斯(Newton-Cotes)法等都是经常采用的方法。它们的基本思想都是将整个积分区间[a,b]分成n个子区间[xi,xi+1]，i=1,2,…,n，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。 2.数值积分的实现方法 基于变步长辛普生法，MATLAB给出了quad函数来求定积分。该函数的调用格式为： [I,n]=quad('fname',a,b,tol,trace) 基于变步长、牛顿－柯特斯(Newton-Cotes)法，MATLAB给出了quadl函数来求定积分。该函数的调用格式为： [I,n]=quadl('fname',a,b,tol,trace) 其中fname是被积函数名。a和b分别是定积分的下限和上限。tol用来控制积分精度，缺省时取tol=0.001。trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，缺省时取trace=0。返回参数I即定积分值，n为被积函数的调用次数。 例： 求函数'exp(-x*x)的定积分，积分下限为0，积分上限为1。 >>fun=inline('exp(-x.*x)','x'); %用内联函数定义被积函数fname 东北大学秦皇岛分校 数值计算课程设计报告 数值积分算法及MATLAB实现 学院数学与统计学院 专业信息与计算科学 学号5133201 姓名陈悦 指导教师姜玉山张建波 成绩 教师评语： 指导教师签字： 2015年07月14日 1 绪论 数值分析是计算数学的一个主要部分，计算数学是数学科学的一个分支，它研究用计算机求解各种数学问题的数值检索方其理论与软件的实现.而数值分析主要研究数值计算. 现科学技术的发展与进步提出了越来越多的复杂的数值计算问题，这些问题的圆满解决已远人工手算所能胜任，必须依靠电子计算机快速准确的数据处理能力.这种用计算机处理数值问题的方法，成为科学计算.今天，科学计算的应用范围非常广泛，天气预报、工程设计、流体计算、经济规划和预测以及国防尖端的一些科研项目，如核武器的研制、导弹和火箭的发射等，始终是科学计算最为活跃的领域. 1.1 数值积分介绍 数值积分是数值分析的重要环节，实际问题当中常常需要计算积分，有些数值方法，如微分方程和积分方程的求解，也都和积分计算相联系. 求某函数的定积分时，在多数情况下，被积函数的原函数很难用初等函数表达出来，因此能够借助微积分学的牛顿-莱布尼兹公式计算定积分的机会是不多的.另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解.由于以上原因，数值积分的理论与方法一直是计算数学研究的基本课题.对微积分学做出杰出贡献的数学大师，如I.牛顿、L.欧拉、C.F.高斯、拉格朗日等人都在数值积分这个领域作出了各自的贡献，并奠定了这个分支的理论基础. 构造数值积分公式最通常的方法是用积分区间上的n 次插值多项式代替被积函数，由此导出的求积公式称为插值型求积公式.特别在节点分布等距的情形称为牛顿-科特斯公式，例如梯形公式(Trapezoidal Approximations)与抛物线公式(Approximations Using Parabolas)就是最基本的近似公式.但它们的精度较差.龙贝格算法是在区间逐次分半过程中，对梯形公式的近似值进行加权平均获得准确程度较高的积分近似值的一种方法，它具有公式简练、计算结果准确、使用方便、稳定性好等优点，因此在等距情形宜采用龙贝格求积公式(Rhomberg Integration).当用不等距节点进行计算时，常用高斯型求积公式计算，它在节点数目相同情况下，准确程度较高，稳定性好，而且还可以计算无穷积分.数值积分还是微分方程数值解法的重要依据.许多重要公式都可以用数值积分方程导出.现探讨数值积分算法以及运用MATLAB软件的具体实现 4.1数值微积分 4.1.1近似数值极限及导数 Matlab 数值计算中，没有求极限指令，也没有求导指令，而是利用差分指令： 用一个简单矩阵表现diff和gradient指令计算方式。 差分： Dx=diff(X) 对向量： Dx=X(2:n)-X(1:n-1) 对矩阵： DX=X(2:n,:)-X(1:n-1,:) 长度小1. DIFF(X), for a vector X, is [X(2)-X(1) X(3)-X(2) ... X(n)-X(n-1)]. DIFF(X), for a matrix X, is the matrix of row differences, （结果缺少一行） [X(2:n,:) - X(1:n-1,:)]. DIFF(X,N,DIM) is the Nth difference function along dimension DIM. If N >= size(X,DIM), DIFF returns an empty array （N阶差分） 梯度： FX=gradient(F) Fx(1)=Fx(2)-Fx(1); F=[1,2,3;4,5,6;7,8,9] Dx=diff(F) （按行） Dx_2=diff(F,1,2) （按列） [FX,FY]=gradient(F) Fx(1)=Fx(2)-Fx(1), Fx(end)=F(end)-F(end-1) FX与F维数相同。 [FX_2,FY_2]=gradient(F,0.5) %采样间隔0.5 即： Fx(1)=(Fx(2)-Fx(1))/2 F = 1 2 3 4 5 6 7 8 9 Dx = 3 3 3 3 3 3 Dx_2 = 1 1 1 1 1 1 FX = 1 1 1 实验一Matlab基本操作与微积分计算 实验目的 1.进一步理解导数概念及其几何意义. 2.学习matlab的求导命令与求导法. 3.通过本实验加深理解积分理论中分割、近似、求和、取极限的思想方法. 4.学习并掌握用matlab求不定积分、定积分、二重积分、曲线积分的方法. 5.学习matlab命令sum、symsum与int. 实验内容 一、变量 1、变量 MA TLAB中变量的命名规则是: （1）变量名必须是不含空格的单个词; （2）变量名区分大小写; （3）变量名最多不超过19个字符; （4）变量名必须以字母打头,之后可以是任意字母、数字或下划线,变量名中不允许使用标点符号. 1、创建简单的数组 x=[a b c d e f ]创建包含指定元素的行向量 x=first:step: last创建从first起,逐步加step计数,last结束的行向量, step缺省默认值为1 x=linspace(first,last,n）创建从first开始,到last结束,有n个元素的行向量 x=logspace(first,last,n）创建从first开始,到last结束,有n个元素的对数分隔行向量. 注:以空格或逗号分隔的元素指定的是不同列的元素,而以分号分隔的元素指定了不同行的元素. 2、数组元素的访问 （1）访问一个元素: x(i)表示访问数组x的第i个元素. （2）访问一块元素: x(a :b :c)表示访问数组x的从第a个元素开始,以步长为b到第c个元素（但 不超过c）,b可以为负数,b缺损时为1. （3）直接使用元素编址序号: x ([a b c d]) 表示提取数组x的第a、b、c、d个元素构成一个新的数组[x (a) x (b) x(c) x(d)]. 3、数组的运算 （1）标量-数组运算 数组对标量的加、减、乘、除、乘方是数组的每个元素对该标量施加相应的加、减、乘、除、乘方运算. 设:a=[a1,a2,…,an], c=标量, 则: a+c=[a1+c,a2+c,…,an+c] a .*c=[a1*c,a2*c,…,an*c] a ./c= [a1/c,a2/c,…,an/c](右除） a .\c= [c/a1,c/a2,…,c/an] (左除） a .^c= [a1^c,a2^c,…,an^c] c .^a= [c^a1,c^a2,…,c^an] （2）数组-数组运算 当两个数组有相同维数时,加、减、乘、除、幂运算可按元素对元素方式进行的,不同大小或维数的数组是不能进行运算的. 设:a=[a1,a2,…,an], b=[b1,b2,…,bn], 则: a +b= [a1+b1,a2+b2,…,an+bn] a .*b= [a1*b1,a2*b2,…,an*bn] a ./b= [a1/b1,a2/b2,…,an/bn] a .\b=[b1/a1,b2/a2,…,bn/an] a .^b=[a1^b1,a2^b2,…,an^bn] 三、矩阵 1、矩阵的建立 矩阵直接输入:从“[ ” 开始,元素之间用逗号“,”(或空格),行之间用分号“;”(或回车),用“ ]”结束. 特殊矩阵的建立: a=[ ] 产生一个空矩阵,当对一项操作无结果时,返回空矩阵,空矩阵的大小为零. b=zeros (m,n) 产生一个m行、n列的零矩阵 c=ones (m,n) 产生一个m行、n列的元素全为1的矩阵 d=eye (m,n) 产生一个m行、n列的单位矩阵 eye (n) %生成n维的单位向量 eye (size (A)) %生成与A同维的单位阵 2、矩阵中元素的操作 （1）矩阵A的第r行A（r,:） （2）矩阵A的第r列A（:,r） （3）依次提取矩阵A的每一列,将A拉伸为一个列向量A（:） （4）取矩阵A的第i1~i2行、第j1~j2列构成新矩阵:A(i1:i2, j1:j2) （5）以逆序提取矩阵A的第i1~i2行,构成新矩阵:A(i2:-1:i1,:） （6）以逆序提取矩阵A的第j1~j2列,构成新矩阵:A(:, j2:-1:j1 ） （7）删除A的第i1~i2行,构成新矩阵:A(i1:i2,:)=[ ] （8）删除A的第j1~j2列,构成新矩阵:A(:, j1:j2)=[ ] （9）将矩阵A和B拼接成新矩阵:[A B];[A;B] 3、矩阵的运算 （1）标量-矩阵运算同标量-数组运算. （2）矩阵-矩阵运算 a. 元素对元素的运算,同数组-数组运算.(A/B %A右除B; B\A%A左除B) b. 矩阵运算: 矩阵加法:A+B 矩阵乘法:A*B 方阵的行列式:det（A） 方阵的逆:inv（A） 用MATLAB 计算多元函数的积分 三重积分的计算最终是化成累次积分来完成的，因此只要能正确的得出各累次积分的积分限，便可在MA TLAB 中通过多次使用int 命令来求得计算结果。但三重积分的积分域Ω是一个三维空间区域，当其形状较复杂时，要确定各累次积分的积分限会遇到一定困难，此时，可以借助MATLAB 的三维绘图命令，先在屏幕上绘出Ω的三维立体图，然后执行命令 rotate3d on ↙ 便可拖动鼠标使Ω的图形在屏幕上作任意的三维旋转，并且可用下述命令将Ω的图形向三个坐标平面进行投影： view(0,0),向XOZ 平面投影； view(90,0),向YOZ 平面投影； view(0,90),向XOY 平面投影. 综合运用上述方法，一般应能正确得出各累次积分的积分限。 例11.6.1计算zdv Ω ???，其中Ω是由圆锥曲面222z x y =+与平面z=1围成的闭区域 解 首先用MA TLAB 来绘制Ω的三维图形，画圆锥曲面的命令可以是： syms x y z ↙ z=sqrt(x^2+y^2); ↙ ezsurf(z,[-1.5,1.5]) ↙ 画第二个曲面之前，为保持先画的图形不会被清除，需要执行命令 hold on ↙ 然后用下述命令就可以将平面z=1与圆锥面的图形画在一个图形窗口内： [x1,y1]=meshgrid(-1.5:1/4:1.5); ↙ z1=ones(size(x1)); ↙ surf(x1,y1,z1) ↙ 于是得到Ω的三维图形如图: 由该图很容易将原三重积分化成累次积分： 111zdv dy -Ω=???? 于是可用下述命令求解此三重积分： clear all ↙ syms x y z ↙ f=z; ↙ f1=int(f,z.,sqrt(x^2+ y^2),1); ↙ f2=int(f1,x,-sqrt(1- y^2), sqrt(1- y^2)); ↙ int(f2,y,-1,1) ↙ ans= 1/4*pi 计算结果为4 π 对于第一类曲线积分和第一类曲面积分，其计算都归结为求解特定形式的定积分和二重积分，因此可完全类似的使用int 命令进行计算，并可用diff 命令求解中间所需的各偏导数。 例11.6.2用MATLAB 求解教材例11.3.1 解 求解过程如下 syms a b t ↙ x=a*cos(t); ↙ y=a*sin(t); ↙ z=b*t; ↙ f=x^2 +y^2+z^2; ↙ xt=diff(x,t); ↙ yt=diff(y,t); ↙ zt=diff(z,t); ↙ int(f*sqrt(xt^2 +yt^2+zt^2),t,0,2*pi) ↙ ans= 2/3*( a^2 +b^2)^1/2*a^2*pi+8/3*( a^2 +b^2)^1/2*b^2*pi^3 对此结果可用factor 命令进行合并化简： factor （ans ） ans= 2/3*( a^2 +b^2)^1/2*pi*(3* a^2 +4*b^2*pi^2) 例11.6.3用MATLAB 求解教材例11.4.1 解 求解过程如下 syms x y z1 z2↙ f= x^2 +y^2; ↙ z1=sqrt(x^2 +y^2); ↙ z2=1; ↙ z1x=diff(z1,x); ↙ z1y=diff(z1,y); ↙ z2x=diff(z2,x); ↙ z2y=diff(z2,y); ↙ M a t l a b数值积分与数值微分 Matlab数值积分 1.一重数值积分的实现方法 变步长辛普森法、高斯-克朗罗德法、梯形积分法 1.1变步长辛普森法 Matlab提供了quad函数和quadl函数用于实现变步长 辛普森法求数值积分.调用格式为: [I,n]=Quad(@fname,a,b,tol,trace) [I,n]=Quadl(@fname,a,b,tol,trace) Fname是函数文件名，a,b分别为积分下限、积分上限； tol为精度控制，默认为1.0×10-6，trace控制是否展 开积分过程，若为0则不展开，非0则展开，默认不展开. 返回值I为积分数值；n为调用函数的次数. --------------------------------------------------------------------- 例如：求 ∫e0.5x sin(x+π )dx 3π 的值. 先建立函数文件 fesin.m function f=fesin(x) f=exp(-0.5*x).*sin(x+(pi/6));再调用quad函数 [I,n]=quad(@fesin,0,3*pi,1e-10) I= 0.9008 n= 365 --------------------------------------------------------------------- 例如：分别用quad函数和quadl函数求积分 ∫e0.5x sin(x+π 6 )dx 3π 的近似值，比较函数调用的次数. 先建立函数文件 fesin.m function f=fesin(x) f=exp(-0.5*x).*sin(x+(pi/6)); formatlong [I,n]=quadl(@fesin,0,3*pi,1e-10) I= n= 198 [I,n]=quad(@fesin,0,3*pi,1e-10) I= n= 365 --------------------------------------------------------------------- 可以发现quadl函数调用原函数的次数比quad少，并 且比quad函数求得的数值解更精确. 1.2高斯-克朗罗德法 matlab 微积分基本运算 §1 解方程和方程组解 1. 线性方程组求解 对于方程 AX = B ,其中 A 是( m ×n )的矩阵有三种情形： 1)当n=m 且A 非奇异时,此方程为“恰定”方程组。 2)当 n > m 时,此方程为“超定”方程组。 3)当n 0.3188 两种方法所求方程组的解相同。 (2)MATLAB 解超定方程AX=B 的方法 对于方程 AX = B ,其中 A 是( m ×n )的矩阵， n > m ，如果A 列满秩，则此方程是没有精确解的。然而在实际工程应用中，求得其最小二乘解也是有意义的。基本解法有： 1)采用求伪逆运算解方程 x=pinv(A)*B 说明：此解为最小二乘解x=inv(A ’*A)*A*B,这里pinv(A) =inv(A ’*A)*A. 2)采用左除运算解方程 x=A\B 例2 “求伪逆”法和“左除”法求下列方程组的解 ??? ??=+=+=+1 221421 221 2121x x x x x x 命令如下： >> a=[1 2;2 4;2 2]; >> b=[1,1,1]'; >> xc=a\b %用左除运算解方程 运行得结果: xc = 0.4000 0.1000 >> xd=pinv(a)*b %用求伪逆运算解方程 运行得结果: xd = 0.4000 0.1000 >> a*xc-b %xc 是否满足方程ax=b 运行得结果: ans = -0.4000 0.2000 0.0000 可见xc 并不是方程的精确解。 (3) MATLAB 解欠定方程AX=B 的方法 欠定方程从理论上说是有无穷多个解的，如果利用求“伪逆”法和“左除”法来求解，只能得到其中一个解。基本方法： 1)采用求伪逆运算解方程 x=pinv(A)*B 2)采用左除运算解方程 1.1 矩阵的表示 1.2 矩阵运算 1.2.14 特殊运算 1．矩阵对角线元素的抽取 函数diag 格式X = diag(v,k) %以向量v的元素作为矩阵X的第k条对角线元素，当k=0时，v为X的主对角线；当k>0时，v为上方第k条对角线；当k<0时，v为下方第k条对角线。 X = diag(v) %以v为主对角线元素，其余元素为0构成X。 v = diag(X,k) %抽取X的第k条对角线元素构成向量v。k=0：抽取主对角线元素；k>0：抽取上方第k条对角线元素；k<0抽取下方第k条对角线元素。 v = diag(X) %抽取主对角线元素构成向量v。 2．上三角阵和下三角阵的抽取 函数tril %取下三角部分 格式L = tril(X) %抽取X的主对角线的下三角部分构成矩阵L L = tril(X,k) %抽取X的第k条对角线的下三角部分；k=0为主对角线；k>0为主对角线以上；k<0为主对角线以下。函数triu %取上三角部分 格式U = triu(X) %抽取X的主对角线的上三角部分构成矩阵U U = triu(X,k) %抽取X的第k条对角线的上三角部分；k=0为主对角线；k>0为主对角线以上；k<0为主对角线以下。3．矩阵的变维 矩阵的变维有两种方法，即用“：”和函数“reshape”，前者主要针对2个已知维数矩阵之间的变维操作；而后者是对于一个矩阵的操作。 （1）“：”变维 （2）Reshape函数变维 格式 B = reshape(A,m,n) %返回以矩阵A的元素构成的m×n矩阵B B = reshape(A,m,n,p,…) %将矩阵A变维为m×n×p×… B = reshape(A,[m n p…]) %同上 B = reshape(A,siz) %由siz决定变维的大小，元素个数与A中元素个数 相同。 （5）复制和平铺矩阵 函数repmat 格式 B = repmat(A,m,n) %将矩阵A复制m×n块，即B由m×n块A平铺而成。 B = repmat(A,[m n]) %与上面一致 B = repmat(A,[m n p…]) %B由m×n×p×…个A块平铺而成 repmat(A,m,n) %当A是一个数a时，该命令产生一个全由a组成的m×n矩阵。 1.3 矩阵分解 1.3.1 Cholesky分解 函数chol 格式R = chol(X) %如果X为n阶对称正定矩阵，则存在一个实的非奇异上三角阵R，满足R'*R = X；若X非正定，则产生错误信息。 [R,p] = chol(X) %不产生任何错误信息，若X为正定阵，则p=0，R与上相同；若X非正定，则p为正整数，R是有序的上三角阵。 1.3.2 LU分解 实验10 数值积分 实验目的： 1．了解数值积分的基本原理； 2．熟练掌握数值积分的MATLAB 实现； 3．会用数值积分方法解决一些实际问题。 实验内容： 积分是数学中的一个基本概念，在实际问题中也有很广泛的应用。同微分一样，在《微积分》中，它也是通过极限定义的，由于实际问题中遇到的函数一般都以列表形式给出，所以常常不能用来直接进行积分。此外有些函数虽然有解析式，但其原函数不是初等函数，所以仍然得不到积分的精确值，如不定积分?1 0 d sin x x x 。这时我们一般考虑用数值方法计算其 近似值，称为数值积分。 10.1 数值微分简介 设函数()y f x =在* x 可导，则其导数为 h x f h x f x f h ) ()(lim )(**0* -+='→ （10.1） 如果函数()y f x =以列表形式给出（见表10-1），则其精确值无法求得，但可由下式求得其近似值 h x f h x f x f ) ()()(*** -+≈' （10.2） 表 10-1 一般的，步长h 越小，所得结果越精确。（10.2）式右端项的分子称为函数()y f x =在 *x 的差分，分母称为自变量在*x 的差分，所以右端项又称为差商。数值微分即用差商近似 代替微商。常用的差商公式为： 000()() ()2f x h f x h f x h +--'≈ （10.3） h y y y x f 243)(2 100-+-≈ ' （10.4） h y y y x f n n n n 234)(12+-≈ '-- （10.5） 其误差均为2 ()O h ，称为统称三点公式。 10.2 数值微分的MATLAB 实现 MATLAB 提供了一个指令求解一阶向前差分，其使用格式为： dx=diff(x) 其中x 是n 维数组，dx 为1n -维数组[]21321,, ,n x x x x x x ---,这样基于两点的数值导 数可通过指令diff(x)/h 实现。对于三点公式，读者可参考例1的M 函数文件diff3.m 。 例1 用三点公式计算()y f x =在=x 1.0,1.2,1.4处的导数值，()f x 的值由下表给 解：建立三点公式的M 函数文件diff3.m 如下： function f=diff3(x,y) n=length(x);h=x(2)-x(1); f(1)=(-3*y(1)+4*y(2)-y(3))/(2*h); for j=2:n-1 f(j)=(y(j+1)-y(j-1))/(2*h); end f(n)=(y(n-2)-4*y(n-1)+3*y(n))/(2*h); 在MATLAB 指令窗中输入指令： x=[1.0,1.1,1.2,1.3,1.4];y=[0.2500,0.2268,0.2066,0.1890,0.1736];diff3(x,y) 运行得各点的导数值为：-0.2470，-0.2170，-0.1890，-0.1650，-0.0014。所以()y f x =在=x 1.0,1.2,1.4处的导数值分别为-0.2470，-0.1890和-0.0014。 对于高阶导数，MATLAB 提供了几个指令借助于样条函数进行求导，详细使用步骤如下： step1：对给定数据点（x,y ），利用指令pp=spline(x,y)，获得三次样条函数数据pp ，供后面ppval 等指令使用。其中，pp 是一个分段多项式所对应的行向量，它包含此多项式的阶数、段数、节点的横坐标值和各段多项式的系数。 step2：对于上面所求的数据向量pp ，利用指令[breaks,coefs,m,n]=unmkpp(pp)进行处理，生成几个有序的分段多项式pp 。 step3：对各个分段多项式pp 的系数，利用函数ppval 生成其相应导数分段多项式的系数，再利用指令mkpp 生成相应的导数分段多项式 step4：将待求点xx 代入此导数多项式，即得样条导数值。 上述过程可建立M 函数文件ppd.m 实现如下： function dy=ppd(pp) [breaks,coefs,m]=unmkpp(pp); 1-1、基本运算与函数 在MATLAB下进行基本数学运算，只需将运算式直接打入提示号（>>）之後，并按入Enter 键即可。例如： >> (5*2+1.3-0.8)*10/25 ans =4.2000 MATLAB会将运算结果直接存入一变数ans，代表MATLAB运算後的答案（Answer）并显示其数值於萤幕上。 小提示： ">>"是MATLAB的提示符号（Prompt），但在PC中文视窗系统下，由於编码方式不同，此提示符号常会消失不见，但这并不会影响到MATLAB的运算结果。 我们也可将上述运算式的结果设定给另一个变数x： x = (5*2+1.3-0.8)*10^2/25 x = 42 此时MATLAB会直接显示x的值。由上例可知，MATLAB认识所有一般常用到的加（+）、减（-）、乘（*）、除（/）的数学运算符号，以及幂次运算（^）。 小提示： MATLAB将所有变数均存成double的形式，所以不需经过变数宣告（Variable declaration）。MATLAB同时也会自动进行记忆体的使用和回收，而不必像C语言,必须由使用者一一指定.这些功能使的MATLAB易学易用，使用者可专心致力於撰写程式，而不必被软体枝节问题所干扰。 若不想让MATLAB每次都显示运算结果，只需在运算式最後加上分号（；）即可，如下例：y = sin(10)*exp(-0.3*4^2); 若要显示变数y的值，直接键入y即可： >>y y =-0.0045 在上例中，sin是正弦函数，exp是指数函数，这些都是MATLAB常用到的数学函数。 下表即为MATLAB常用的基本数学函数及三角函数： 小整理：MATLAB常用的基本数学函数 abs(x)：纯量的绝对值或向量的长度 angle(z)：复数z的相角(Phase angle) sqrt(x)：开平方 real(z)：复数z的实部 imag(z)：复数z的虚部 conj(z)：复数z的共轭复数 round(x)：四舍五入至最近整数 fix(x)：无论正负，舍去小数至最近整数 floor(x)：地板函数，即舍去正小数至最近整数 ceil(x)：天花板函数，即加入正小数至最近整数 rat(x)：将实数x化为分数表示 rats(x)：将实数x化为多项分数展开 sign(x)：符号函数 (Signum function)。 当x<0时，sign(x)=-1； 当x=0时，sign(x)=0; 当x>0时，sign(x)=1。 > 小整理：MATLAB常用的三角函数 sin(x)：正弦函数 科学计算—理论、方法 及其基于MATLAB 的实现与分析 数值微积分 §1 数值微分 对于给定的函数()x f y =，如果 １、()x f y =的函数关系式比较复杂时； ２、 ()x f y =未知，而仅仅知道()x f y =在1+n 个相异点k x ， n k ,,1,0Λ=处 的函数值k y ； 则希望能用相对简单的计算方法，求得()x f y =导数的（近似）值。 基于上述考虑，选择的方法之一是利用函数()x f y =的插值多项式的导数作为函数()x f y =导数的近似值，例如Lagrange 插值多项式，由于 ()()()()∑==≈n k k k n x f x l x L x f 0 （１） 因而有 ()()x L x f n '≈' （２） 这里需要说明一点的是，尽管()x f 和()x L n 的函数值可能相差不多，但是它们的导数有可能相差很大，如下面的例子 例1： 考虑函数()2 2511x x f += 在区间[-1，1]的插值问题，取 区间[-1，1]的11个点 k x k ?+-=2.01，10,,1,0Λ=k ，作函数 ()2 2511x x f += 的10次插值多项式： ()18552.163597.1234338.3819095.4949417.220246810+-+-+-=x x x x x x L n 函数()x f 和插值多项式()x L n 的导数分别为 ()() 2 225150x x dx x df +-= ()x x x x x dx x dLn 7.334.4936.22883.39594.22093579-+-+-= 对函数()x f 和插值多项式()x L n 及其导数分别比较，结果如图所示： Derivative_Runge 实验一 MATLAB 运算基础 一、实验目的 1.熟悉启动和退出MATLAB 的方法； 2.熟悉MATLAB 命令窗口的组成； 3.掌握建立矩阵的方法； 4.掌握MATLAB 各种表达式的书写规则以及常用函数的使用。 二、实验内容 1.先求下列表达式的值，然后显示MATLAB 工作空间的使用情况并保存全部变量。 ⑴21185sin 2e z +?=； >> z1=2*sin(85*pi/180)/(1+exp(2)) z1 = 0.2375 ⑵)1ln(2122x x z ++=，其中?? ????-+=545.0212i x ； >> x=[2 1+2i;-0.45 5]; >> z2=1/2*log(x+sqrt(1+x^2)) z2 = 0.7114 - 0.0253i 0.8968 + 0.3658i 0.2139 + 0.9343i 1.1541 - 0.0044i ⑶0.3,9.2,8.2,,8.2,9.2,0.3,2 3.0ln )3.0sin(23.03.03 ---=+++-=-a a a e e z a a >> a=(-3.0:0.1:3.0); >> z3=(exp(0.3.*a)-exp(-0.3.*a))./2.*sin(a+0.3)+log((0.3+a)./2) z3 = Columns 1 through 3 0.7388 + 3.1416i 0.7696 + 3.1416i 0.7871 + 3.1416i Columns 4 through 6 0.7913 + 3.1416i 0.7822 + 3.1416i 0.7602 + 3.1416i Columns 7 through 9 MatLab & 数学建模 第一讲简介及基本运算 一、简介 MATLAB名字由MATrix和 LABoratory 两词的前三个字母组合而成。那是20世纪七十年代后期的事：时任美国新墨西哥大学计算机科学系主任的Cleve Moler教授出于减轻学生编程负担的动机，为学生设计了一组调用LINPACK和EISPACK库程序的“通俗易用”的接口，此即用FORTRAN编写的萌芽状态的MATLAB。 经几年的校际流传，在Little的推动下，由Little、Moler、Steve Bangert合作，于1984年成立了MathWorks公司，并把MATLAB正式推向市场。从这时起，MATLAB的内核采用C语言编写，而且除原有的数值计算能力外，还新增了数据图视功能。 MATLAB以商品形式出现后，仅短短几年，就以其良好的开放性和运行的可靠性，使原先控制领域里的封闭式软件包（如英国的UMIST，瑞典的LUND和SIMNON，德国的KEDDC）纷纷淘汰，而改以MATLAB为平台加以重建。在时间进入20世纪九十年代的时候，MATLAB 已经成为国际控制界公认的标准计算软件。 在欧美大学里，诸如应用代数、数理统计、自动控制、数字信号处理、模拟与数字通信、时间序列分析、动态系统仿真等课程的教科书都把MATLAB作为内容。这几乎成了九十年代教科书与旧版书籍的区别性标志。在那里，MATLAB是攻读学位的大学生、硕士生、博士生必须掌握的基本工具。 在国际学术界，MATLAB已经被确认为准确、可靠的科学计算标准软件。在许多国际一流学术刊物上，（尤其是信息科学刊物），都可以看到MATLAB的应用。 在设计研究单位和工业部门，MATLAB被认作进行高效研究、开发的首选软件工具。如美国National Instruments公司信号测量、分析软件LabVIEW，Cadence公司信号和通信分析设计软件SPW等，或者直接建筑在MATLAB之上，或者以MATLAB为主要支撑。又如HP公司的VXI硬件，TM公司的DSP，Gage公司的各种硬卡、仪器等都接受MATLAB的支持。 MATLAB具有用法简易、可灵活运用、程式结构强又兼具延展性。以下为其几个特色： ?功能强的数值运算 - 在MATLAB环境中，有超过500种数学、统计、科学及工程方面的函数可使用，函数的标示自然，使得问题和解答像数学式子一般简单明了，让使用者可全力发挥在解题方面，而非浪费在电脑操作上。 ?先进的资料视觉化功能 - MATLAB的物件导向图形架构让使用者可执行视觉数据分，并制作高品质的图形，完成科学性或工程性图文并茂的文章。 ?高阶但简单的程式环境 - 作为一种直译式的程式语言，MATLAB容许使用者在短时间内写完程式，所花的时间约为用 FORTRAN 或 C 的几分之一，而且不需要编译 (compile)及联结 (link) 即能执行，同时包含了更多及更容易使用的内建功能。 ?开放及可延伸的架构 - MATLAB容许使用者接触它大多数的数学原使码，检视运算法，更改现存函数，甚至加入自己的函数使 MATLAB成为使用者所须要的环境。 ?丰富的程式工具箱 - MATLAB的程式工具箱融合了套装前软体的优点，与一个灵活的开放但容易操作之环境，这些工具箱提供了使用者在特别应用领域所需之许多函数。现有工具箱有：符号运算（利用Maple V的计算核心执行）、影像处理、统计分析、讯号处理、神经网路、模拟分析、控制系统、即时控制、系统确认、强建控制、弧线分析、最佳化、模糊逻辑、mu分析及合成、化学计量分析。 二、MatLab界面MATLAB运算基础第2章答案
用递推公式计算定积分(matlab版)
matlab数值微积分与方程数值求解
实验一 基本信号在MATLAB中的表示和运算
matlab求定积分之实例说明
数值积分算法与MATLAB实现陈悦5133201讲解
MATLAB数值微积分
实验一B Matlab基本操作与微积分计算
用MATLAB算多元函数积分
Matlab数值积分与数值微分
matlab微积分基本运算
matlab中的矩阵的基本运算命令
数值积分的matlab实现
matlab基本运算与函数
第七章 基于MATLAB的科学计算—数值微积分1
实验一 MATLAB运算基础
MatLab基本运算
相关主题
文本预览