15. 数值计算—高数篇

一、求极限

limit(f,x,a)——求极限

limit(f,x,a,'right')——求右极限

limit(f,x,a,'left')——求左极限

例1求

代码：

syms x;

y=(3*x^2+5)/(5*x+3)*sin(2/x);

limit(y,x,inf)

运行结果：ans =6/5

注：Matlab求二元函数的极限，是用嵌套limit函数实现的，相当于求的是累次极限，需要注意：有时候累次极限并不等于极限。

例2求

代码：

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 求的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设，求

代码：

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))*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求函数在区间(-2,4)上的极值

代码：

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.0000min =-47.0000

x1 = -1.0000max = -17.0000

2. 多元函数极值

[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.0000

f1 = 2.7969e-010

四、求不定积分与定积分

1. 符号积分

int(f,x)——求f(x)关于x的不定积分

int(f,x,a,b)——求f(x)关于x的从a到b的定积分

例7求积分和

代码：

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求二重积分, 先化为累次积分：

原式=

代码：

syms x y;

int(int(1+x+y,y,-sqrt(1-x^2),sqrt(1-x^2)),x,-1,1)

运行结果：ans =pi

3. 数值积分

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求

代码：

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求在x=4处展开到2阶泰勒式，在的傅里叶展开。

代码：

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)——

例12求级数

syms n;

symsum((-1)^(n+1)/(n+1)^2,n,1,inf)

运行结果：ans =1 - pi^2/12

七、代数方程

1. 求代数方程的解析解

solve(‘eq1’,’eq2’,…,’var1’,’var2’,…)

例13解方程的x和b，以及方程组

代码：

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/3

2. 非线性方程（组）数值解

fsolve(f,x0)

——求方程f(x)=0在x0附近的近似解，也可以解方程组

注：一元连续函数的根，可以用fzero(f,x0)

例14求解方程。

代码：

f=@(x) x-exp(-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.5671

fval = 1.0e-006 * -0.4059

-0.4059

八、常微分方程（组）

1. 求解析解

dsolve(‘eq1’,’eq2’,…,’cond1’,’cond2’,…,’t’)

默认自变量为t，cond1,2…为初值条件，若有足够初值条件，则得到特解；否则得到通解。若解不出解析解，只能用ode23或ode45求数值解。用Dy, D2y,…表示；用D2y(e)=a表示.

例16求解微分方程

代码：

y1=dsolve('y*D2y-Dy^2=0','x')

运行结果：y1 = C4

exp(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. 求数值集（利用求解器）

实际问题中，许多常微分方程（组）是求不出解析解的，Matlab提供了多个求数值集的求解器solver.

调用格式：

[T,X]=solver(odefun, 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)

Matlab笔记 - 数值计算—高数篇015