(完整版)Matlab学习系列15.数值计算—高数篇

15. 数值计算—高数篇一、求极限limit(f,x,a)——求极限 lim ()x af x →limit(f,x,a,'right')——求右极限 lim ()x af x +→limit(f,x,a,'left')——求左极限 lim ()x af x -→例1 求 2352lim sin 53x x x x→∞++代码：syms x;y=(3*x^2+5)/(5*x+3)*sin(2/x); limit(y,x,inf)运行结果：ans = 6/5注：Matlab 求二元函数的极限，是用嵌套limit 函数实现的，相当于求的是累次极限，需要注意：有时候累次极限并不等于极限。

例2 求 30lim 2x x xx a b +→⎛⎫+ ⎪⎝⎭代码：syms x a b;y=((a^x+b^x)/2)^(3/x); limit(y,x,0,'right')运行结果：ans = a^(3/2)*b^(3/2)二、求导diff(f,x,n)——求函数f 关于x 的n 阶导数，默认n=1例3 求1sin 1cos xy x+=+的1阶导数，并绘图代码：syms x a b;y=((a^x+b^x)/2)^(3/x); limit(y,x,0,'right')运行结果：y1 = cos(x)/(cos(x) + 1) + (sin(x)*(sin(x) + 1))/(cos(x) + 1)^2例4 设sin xyz e=，求2,,z z z x y x y∂∂∂∂∂∂∂代码：syms x y;z=exp(sin(x*y)); zx=diff(z,x) zy=diff(z,y)zxy=diff(zx,y) % 也等于diff(zy,x)-55510152025x(sin(x) + 1)/(cos(x) + 1)-505-20-15-10-50510152025xcos(x)/(cos(x) + 1) + (sin(x) (sin(x) + 1))/(cos(x) + 1)2运行结果：zx = y*exp(sin(x*y))*cos(x*y)zy = x*exp(sin(x*y))*cos(x*y)zxy = x*y*exp(sin(x*y))*cos(x*y)^2 + exp(sin(x*y))*cos(x*y)–x*y*exp(sin(x*y))*sin(x*y)三、求极值1. 一元函数极值[x0,min]=fminbnd(f, a, b)——返回函数f(x)在区间(a,b)上的极小值点和极小值例5求函数32=--+在区间(-2,4)上的极值()26187f x x x x代码：f=@(x) 2*x^3-6*x^2-18*x+7;g=@(x) -2*x^3+6*x^2+18*x-7;[x0,min]=fminbnd(f,-2,4)[x1,max]=fminbnd(g,-2,4)fplot(f,[-2,4]);运行结果：x0 = 3.0000 min = -47.0000x1 = -1.0000 max = -17.00002. 多元函数极值[X1,f1]=fminunc(f,X0)——处理连续情形[X1,f1]=fminsearch(f,X0)——可以处理不连续情形二者用法相同，返回极小值点和极小值，其中X 为初始点。

例6 求222(,)(1)100()f x y x y x =-+-的极小值代码：f=@(x) (1-x(1))^2+100*(x(2)-x(1)^2)^2; x0=[-5 -2];[x1,f1]=fminsearch(f,x0)运行结果：x1 = 1.0000 1.0000f1 = 2.7969e -010四、求不定积分与定积分1. 符号积分int(f,x)——求f(x)关于x 的不定积分01020int(f,x,a,b)——求f(x)关于x 的从a 到b 的定积分 例7 求积分2ln d x ax x -⎰和21ln d x a x x+∞-⎰ 代码：syms x a;int((log(x)-a)/x^2,x)int((log(x)-a)/x^2,x,1,inf)运行结果：ans = -(log(x) - a + 1)/x ans = 1 – a注：不定积分的结果是忽略任意常数C 的。

2. 二重积分可以化为累次积分，再用两次int 函数实现。

例8 求二重积分221(1)d d x y x y x y +≤++⎰⎰, 先化为累次积分：原式=11d )d x x y y -++⎰代码：syms x y;int(int(1+x+y,y,-sqrt(1-x^2),sqrt(1-x^2)),x,-1,1)运行结果：ans = pi3. 数值积分quad(f, a, b)——辛普森法定积分，默认误差为10-6，低精度的非光滑曲线计算中是最有效；quadl(f, a, b)——Lobatto 法定积分，在高精度的光滑曲线计算中更为高效；quad2d(f, a, b, c, d)——二重积分，其中f(x,y)为二元函数，[a, b]为x 的范围，[c(x), d(x)]为y 的范围； 例9 求1021ln 1d x x x-⎰代码：f=@(x) (log(x)-1)./x.^2; % 注意’.’不能忽略 y=quad(f,1,10)运行结果：y = -0.2303例10 用数值积分法求解例8.代码：quad2d(@(x,y) 1+x+y, -1, 1, @(x) -sqrt(1-x.^2), @(x) sqrt(1-x.^2), 'AbsTol', 1e-12) % 注意点运算运行结果：ans = 3.1416或者用两次quad 函数，中间需要用arrayfun 函数做向量值化处理，该方法可以推广到三重积分。

quad(@(x) arrayfun(@(xx) quad(@(y) 1+xx+y,-sqrt(1-xx.^2),sqrt(1-xx.^2)),x), -1,1)程序说明：先对f(x,y)关于y 从-sqrt(1-x.^2)到sqrt(1-x.^2)]做一次积分，为了后面使用变量名x ，这里先用xx ，得到一个关于xx 的函数（只能接受自变量为标量xx ）：quad(@(y) 1+xx+y, -sqrt(1-xx.^2),sqrt(1-xx.^2))然后用arrayfun 函数把上一步得到的xx 的函数，处理成能接受向量值x （是个x 的函数）：@(x)arrayfun(@(xx)quad(@(y) 1+xx+y, -sqrt(1-xx.^2),sqrt(1-xx.^2)), x)最后，再关于x 做一次积分。

五、泰勒级数、傅里叶级数展开taylor(f,n,x,a)——将函数f(x)在x=a 点处展开为n -1阶泰勒级数 fseries(f,x,n,a,b)——将函数f(x)在区间(a,b)展开n 项傅里叶级数 注：Matlab 未提供傅里叶级数展开函数，fseries 函数来自论坛。

例11 求()10af x x =-在x=4处展开到2阶泰勒式，2()g x x x =+在[,]ππ-的傅里叶展开。

代码：syms a x; f=a/(x-10);y1=taylor(f,3,x,4) g=x^2+x;[an,bn,f]=fseries(g,x,3,-pi,pi)运行结果：y1 = - a/6 - (a*(x - 4))/36 - (a*(x - 4)^2)/216 an = [ (2*pi^2)/3, -4, 1, -4/9] bn = [ 2, -1, 2/3]f = cos(2*x) - (4*cos(3*x))/9 - sin(2*x) + (2*sin(3*x))/3 - 4*cos(x) + 2*sin(x) + pi^2/3六、求级数symsum(f, k, m,n)——()nk mf k =∑例12 求级数121(1)(1)n n n +∞=-+∑ syms n;symsum((-1)^(n+1)/(n+1)^2,n,1,inf)运行结果：ans = 1 - pi^2/12七、代数方程1. 求代数方程的解析解solve(‘eq1’,’eq2’,…,’var1’,’var2’,…)例13 解方程20ax bx c ++=的x 和b ，以及方程组1115x y x y +=⎧⎨-=⎩代码：syms a b c x;solve('a*x^2+b*x+c','x') solve('a*x^2+b*x+c','b')[x,y]=solve('x+y=1','x-11*y=5','x','y')运行结果：ans = -(b + (b^2 - 4*a*c)^(1/2))/(2*a) -(b - (b^2 - 4*a*c)^(1/2))/(2*a) ans = -(a*x^2 + c)/x x = 4/3 y = -1/32. 非线性方程（组）数值解fsolve(f,x0)——求方程f(x)=0在x0附近的近似解，也可以解方程组 注：一元连续函数的根，可以用fzero(f,x0)例14 求解方程0x x e --=。

代码：f=@(x) x-exp(-x); x1=fsolve(f,0)运行结果：x1 = 0.5671例15 求解方程组 2020x yx y e x y e --⎧--=⎨-+-=⎩ 代码：F=@(x) [2*x(1)-x(2)-exp(-x(1)); -x(1)+2*x(2)-exp(-x(2))]; [x,fval]=fsolve(F,[-5;-5])运行结果：x = 0.56710.5671fval = 1.0e -006 * -0.4059 -0.4059 八、常微分方程（组）1. 求解析解dsolve(‘eq1’,’eq2’,…,’cond1’,’cond2’,…,’t’)默认自变量为t ，cond1,2…为初值条件，若有足够初值条件，则得到特解；否则得到通解。

若解不出解析解，只能用ode23或ode45求数值解。

用Dy, D2y,…表示','',...y y ；用D2y(e)=a 表示''()|x e y x a ==. 例16 求解微分方程 2''''0y y y -=代码：y1=dsolve('y*D2y-Dy^2=0','x')运行结果：y1 = C4exp(C3 - C2*x) （注：两个解）例17 求解微分方程组 020d d 210cos , |2d d d d 24, |0d d t t t xy x t x t tx y y e y t t=-=⎧--==⎪⎪⎨⎪++==⎪⎩代码：[X,Y]=dsolve('Dx+2*x-Dy=10*cos(t)','Dx+Dy+2*y=4*exp(-2*t)','x(0)=2','y(0)=0')运行结果：X = 4*cos(t) - 2/exp(2*t) + 3*sin(t) - (2*sin(t))/exp(t) Y = sin(t) - 2*cos(t) + (2*cos(t))/exp(t)2. 求数值集（利用求解器）实际问题中，许多常微分方程（组）是求不出解析解的，Matlab 提供了多个求数值集的求解器solver.调用格式：[T,X]=solver(odefun, tspan, X0)其中，tspan 为求解区间；X 0为初值条件向量；先改写高阶微分方程()(,',...,,)0n F y y y t =做变量代换处理：令(1)12,',...,n n y y y y y y -===，得到'11'22'(,)(,)'(,)nn f t Y y f t Y y Y f t Y y ⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪== ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭M M ， 120(1)(0)(0)(0)'(0)(0)(0)n n y y y y Y y y -⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪== ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭M M 例18 求解描述振荡器的经典Ver der Pol 微分方程（取7μ=）：222d d (1)0d d y y y y t tμ--+=， (0)1y =， '(0)0y = 做变量代换处理，12d ,d y x y x t==，则 1222121d d d 7(1)d x x t x x x x t ⎧=⎪⎪⎨⎪=--⎪⎩代码：先编写VDP.m 函数function fy=VDP(t,x)fy=[x(2);7*(1-x(1)^2)*x(2)-x(1)];end主程序：Y0=[1;0];[t,x]=ode45('VDP',[0,40],Y0);y=x(:,1);dy=x(:,2);plot(t,y,t,dy); legend('y','y''');运行结果：注：想得到解y(t)在t0处的值，可以[t,x]=ode45('VDP',[0,t0],Y0);y=x(:,1);y(end) 0510152025303540-15-10-551015。