Matlab笔记数值计算高数篇015.docx

1、Matlab笔记数值计算高数篇01515. 数值计算高数篇一、求极限limit(f,x,a)求极限 limit(f,x,a,right)求右极限limit(f,x,a,left)求左极限例1求代码：syms x;y=(3*x2+5)/(5*x+3)*sin(2/x);limit(y,x,inf)运行结果：ans =6/5注：Matlab求二元函数的极限，是用嵌套limit函数实现的，相当于求的是累次极限，需要注意：有时候累次极限并不等于极限。例2求代码：syms x a b;y=(ax+bx)/2)(3/x);limit(y,x,0,right)运行结果：ans =a(3/2)*b(3/2)二

2、、求导 diff(f,x,n)求函数f关于x的n阶导数，默认n=1例3 求的1阶导数，并绘图代码：syms x a b;y=(ax+bx)/2)(3/x);limit(y,x,0,right)运行结果：y1 =cos(x)/(cos(x) + 1) + (sin(x)*(sin(x) + 1)/(cos(x) + 1)2例4设，求代码：syms x y;z=exp(sin(x*y);zx=diff(z,x)zy=diff(z,y)zxy=diff(zx,y) % 也等于diff(zy,x)运行结果：zx =y*exp(sin(x*y)*cos(x*y)zy =x*exp(sin(x*y)*co

3、s(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求函数在区间(-2,4)上的极值代码：f=(x) 2*x3-6*x2-18*x+7;g=(x) -2*x3+6*x2+18*x-7;x0,min=fminbnd(f,-2,4)x1,max=fminbnd(g,-2,4)fplot(f,-2,4);运行结果：x0 = 3.0000min =-47

4、.0000x1 = -1.0000max = -17.00002. 多元函数极值X1,f1=fminunc(f,X0)处理连续情形X1,f1=fminsearch(f,X0)可以处理不连续情形二者用法相同，返回极小值点和极小值，其中X为初始点。例6 求的极小值代码：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的不定积分int(f,x,a,b)求f(x)关于x的从a到b

5、的定积分例7求积分和代码：syms x a;int(log(x)-a)/x2,x)int(log(x)-a)/x2,x,1,inf)运行结果：ans =-(log(x) - a + 1)/xans =1 a注：不定积分的结果是忽略任意常数C的。2. 二重积分可以化为累次积分，再用两次int函数实现。例8求二重积分, 先化为累次积分：原式=代码：syms x y;int(int(1+x+y,y,-sqrt(1-x2),sqrt(1-x2),x,-1,1)运行结果：ans =pi3. 数值积分quad(f,a,b)辛普森法定积分，默认误差为10-6，低精度的非光滑曲线计算中是最有效；quadl(f

6、, a, b)Lobatto法定积分，在高精度的光滑曲线计算中更为高效；quad2d(f, a, b, c, d)二重积分，其中f(x,y)为二元函数，a, b为x的范围，c(x), d(x)为y的范围；例9求代码：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函数，中

7、间需要用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(x

8、x) 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求在x=4处展开到2阶泰勒式，在的傅里叶展开。代码：syms a x;f=a/(x-10);y1=taylor(f,3,x,4)g=x2+x;an,bn,f=fseries(g,x,3,-pi,pi)运行结果：y

9、1 =- a/6 - (a*(x - 4)/36 - (a*(x - 4)2)/216an = (2*pi2)/3, -4, 1, -4/9bn = 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) + pi2/3六、求级数symsum(f, k, m,n)例12求级数syms n;symsum(-1)(n+1)/(n+1)2,n,1,inf)运行结果：ans =1 - pi2/12七、代数方程1. 求代数方程的解析解solve(eq1,eq2,var1,var2,)例

10、13解方程的x和b，以及方程组代码：syms a b c x;solve(a*x2+b*x+c,x)solve(a*x2+b*x+c,b)x,y=solve(x+y=1,x-11*y=5,x,y)运行结果：ans = -(b + (b2 - 4*a*c)(1/2)/(2*a)-(b - (b2 - 4*a*c)(1/2)/(2*a)ans =-(a*x2 + c)/xx =4/3 y =-1/32. 非线性方程（组）数值解fsolve(f,x0)求方程f(x)=0在x0附近的近似解，也可以解方程组注：一元连续函数的根，可以用fzero(f,x0)例14求解方程。代码：f=(x) x-exp(-

11、x);x1=fsolve(f,0)运行结果：x1 = 0.5671例15求解方程组代码：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.5671 0.5671fval = 1.0e-006 * -0.4059 -0.4059八、常微分方程（组）1. 求解析解dsolve(eq1,eq2,cond1,cond2,t)默认自变量为t，cond1,2为初值条件，若有足够初值条件，则得到特解；否则得到通解。若解不出解析解，只能用ode23或ode45求数值解。用Dy, D2y,表示

12、；用D2y(e)=a表示.例16求解微分方程代码：y1=dsolve(y*D2y-Dy2=0,x)运行结果：y1 = C4exp(C3 - C2*x)（注：两个解）例17求解微分方程组代码：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. 求数值集（利用求解器）实际问题中，许多常微分方程（组）是求不出解析解的，M

13、atlab提供了多个求数值集的求解器solver.求解器ODE类型特点说明ode45非刚性一步算法，4,5阶Runge-Kutta方法累积截断误差大部分场合的首选算法ode23非刚性一步算法，2,3阶Runge-Kutta方法累积截断误差使用于精度较低的情形ode113非刚性多步法，Adams算法，高低精度均可达到计算时间比ode45短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法，Gears反向数值积分，精度中等若ode45失效时，可尝试使用ode23s刚性一步法，2阶Rosebrock算法，低精度。当精度较低时，计算时间比ode15s短调用格式：T,X=solver(o

14、defun, tspan, X0)其中，tspan为求解区间；X0为初值条件向量；先改写高阶微分方程做变量代换处理：令，得到， 例18求解描述振荡器的经典Ver der Pol微分方程（取）：， 做变量代换处理，则代码：先编写VDP.m函数functionfy=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)