Matlab学习系列15数值计算-高数篇

15.数值计算一高数篇

一、求极限

limit(f,x,a)------求极限limf(x)

x-^a

limit(f,x,a,Tigh£)-----求右极限limJ\x)

limit(f,x,a,'left,)-----求左极限lim/(x)

例1求+^sin—

xtb5x+3x

代码：

symsx;

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

limit(y,x,inf)

运行结果：ans=6/5

注：Matlab求二元函数的极限，是用嵌套limit函数实现的，相

当于求的是累次极限，需要注意：有时候累次极限并不等于极限。

mI,+1-(优+"丫

例2求lim---------

52)

代码：

symsxab;

y=((aAx+bAx)/2)A(3/x);

limit(yzx,0,'right*)

运行结果：ans=aA(3/2)*bA(3/2)

二、求导

diff(f,x,n)一一求函数f关于x的n阶导数，默认n=l

例3求歹=上@”的1阶导数，并绘图

1+COSX

代码：

symsxab;

y=((aAx+bAx)/2)A(3/x);

limit(yzx,0,'right*)

运行结果：

yl=cos(x)/(cos(x)+1)+(sin(x)*(sin(x)+1))/(cos(x)+1)A2

(sin(x)+1)/(cos(x)+1)cos(x)/(cos(x)+1)+Sin(x)(sin(x)+1))/(cos(x)+1)*

2525

20

20

15

10

15

5

0

10

-5

5-10

-15

0

-20

-505-505

XX

卡dz8zd2z

例4设Z=求丁,一，----

8xdydxoy

代码:

symsxy;

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)A2+exp(sin(x*y))*cos(x*y)

-x*y*cxp(sin(x*y))*sin(x*y)

三、求极值

1.一元函数极值

[xO,minJ=frninbnd(f,a,b)

----返回函数f(x)在区间(a,b)上的极小值点和极小值

例5求函数/(x)=2/_6/_18x+7在区间(・2,4)上的极值

代码：

f=@(x)2*xA3-6*xA2-18*x+7;

g=0(x)-2*xA3+6*xA2+18*x-7;

[xO,min]=fminbnd(f,-2,4)

[xlzmax]=fminbnd(g,-2,4)

fplojf,[-2,4]);

运行结果：xO=3.0000min=-47.0000

xl=-1.0000max=-17.0000

20

10

0

-10

-20

-30

-40

-50

-2-101234

2.多元函数极值

[Xl,fl]=fminunc(f,X0)——处理连续情形

[X1,f1]=fminsearch(f,XO)-----可以处理不连续情形

二者用法相同，返回极小值点和极小值，其中X为初始点。

例6求/(x,y)=(l-x)2+100Cy--)2的极小值

代码：

f=@(x)(1-x(1))^2+100*(X(2)-X(1)A2)A2;

x0=[-5-2];

[xl,f1]=fminsearch(f,xO)

运行结果：xl=1.00001.0000

fl=2.7969e-010

四、求不定积分与定积分

1.符号积分

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

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

例7求积分］中心和JJ卡dx

代码：

symsxa;

int((log(x)-a)/xA2,x)

int((log(x)-a)/xA2,x,1,inf)

运行结果：ans=-(log(x)-a+l)/x

ans=1-a

注：不定积分的结果是忽略任意常数C的。

2.二重积分

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

例8求二重积分“(l+x+y)&dy,先化为累次积分：

x2+y25l

原式可£(l+x+y)dy

代码：

symsxy;

int(int(1+x+y,y,-sqrt(l-xA2),sqrt(l-xA2)),x,-1,1)

运行结果：ans=pi

3.数值积分

quad(f,a,b)——辛普森法定积分，默认误差为1OR低精度的

非光滑曲线计算中是最有效；

quadl(f,a,b)------Lobatto法定积分，在高精度的光滑曲线计算

中更为局效；

quad2d(f,a,b,c：d)---二重积分，其中f(x,y)为二元函数，［a,b］

为x的范围，［c(x),d(x)］为y的范围；

例9求『写几

代码：

f=0(x)(log(x)-l)./x.A2;茗注意'不能忽略

y=quad(fz1,10)

运行结果：y=-0.2303

例10用数值积分法求解例8.

代码：

quad2d(@(x,y)1+x+y,-1,1,@(x)-sqrt(1-x.A2),0(x)

sqrt(l-x."2)z•Abslol',le-12)之注意点运算

运行结果：ans=3.1416

或者用两次quad函数，中间需要用arrayfun函数做向量值化处

理，该方法可以推广到三重积分。

quad(0(x)arrayfun(0(xx)quad(0(y)1+xx+y,

AA

-sqrt(1-xx.2),sqrt(1-xx.2))zx),-1,1)

程序说明：

先对f(x,y)关于y从-sqrt(l-x.^2)至sqrt(l-x,八2)］做一次积

分，为了后面使用变量名x,这里先用xx,得到一个关于XX的函数

(只能接受自变量为标量xx)：

quad(0(y)1+xx+y,-sqrt(1-xx.A2),sqrt(1-xx.A2))

然后用arrayfun函数把上一步得到的xx的函数，处理成能接受

向量值x(是个x的函数):

@(X)

arrayfun(@(xx)

quad(0(y)1+xx+y,-sqrt(1-xx.A2),sqrt(1-xx.A2)),

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未提供傅里叶级数展开函数，fserics函数来自论坛。

例11求/(、)=―-—在x=4处展开到2阶泰勒式，g(x)=%2+x在

x-10

[-乃的傅里叶展开。

代码：

symsax;

f=a/(x-10);

yl=taylor(f,3,xz4)

g=xA2+x;

[an,bn,f]=fseries(gzx,3,-pi,pi)

运行结果：

yl=-a/6-(a*(x-4))/36-(a*(x-4)A2)/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)+piA2/3

六、求级数

symsum(f)k,m,n)ifw

k-m

8/_i\n+l

例12求级数£匚二

£(〃+l)2

symsn;

symsum((-1)八(r.+l)/(n+1)八2,n,l,inf)

运行结果：ans=1-piA2/12

七、代数方程

1.求代数方程的解析解

solveCeql,,,eq2\...,,varr,,var2,,...)

x+y=\

例13解方程+/zx+c=O的x和b,以及方程组

x-\\y=5

代码：

symsabcx;

solve{'a*xA2+b*x+c*,'x*)

solve(*a*xA2+b*x+c*,'b*)

[x,y]=solve(*x+y=l*,'x-ll*y=5',*x',*y')

运行结果：ans=-(b+(bA2-4*a*c)A(l/2))/(2*a)

-(b-(bA2-4*a*c)A(l/2))/(2*a)

ans=-(a*xA2+c)/x

x=4/3y=-l/3

2.非线性方程(组)数值解

fsolve(f,xO)

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

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

例14求解方程x-"X=o。

代码：

f=0(x)x-exp(-x);

xl=fsolve(f,0)

运行结果：xl=0.5671

例15求解方程组[2%一尸。一、二°

-x+2y-e-=0

代码：

F=0(x)[2*x(l)-x(2)-exp(-x(l));

-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(4eq1\,eq2,,.../condl\,cond2\...,,t,)

默认自变量为t,condl,2...为初值条件，若有足够初值条件，则

得到特解；否则得到通解。若解不出解析解，只能用ode23或。de45

求数值解。用Dy,D2y,...表示…;用D2y(e尸a表示y"(x)二=衣

例16求解微分方程-尸2=0

代码：

yl=dsolve(*y*D2y-Dy^2=0*,'x1)

运行结果：yl=C4

exp(C3-C2*x)(注：两个解)

drcdy

-----2x-----=lOcos/,了l-o=2

dtdt

例17求解微分方程组

dvdy2f

—+—+2y=4e~,“0=

&dz。

代码：

[X,Y]=dsolve(1Dx+2*x-Dy=10*cos(t)',1Dx+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.

求解器ODE类型特点说明

•步算法，4,5阶Runge-Kutta

ode45非刚性大部分场合的首选算法

方法累积截断误差(右)3

一步算法,2,3阶Runge-Kutta

ode23非刚性使用于精度较低的情形

方法累积截断误差(以)3

多步法，Adams算法,高低精

odel13非刚性计算时间比ode45短

度均可达到1°-3〜1°^

ode23t适度刚性采用梯形算法适度刚性情形

多步法，Gear's反向若ode45失效时，

ode15s刚性

数值积分，精度中等可尝试使用

ode23s刚性一步法，2阶Rosebrock算法，当精度较低时，

低精度。计算时间比ode15s短

调用格式：

[T,X]=solver(odefun,tspan,XO)

其中，tspan为求解区间；X。为初值条件向量；先改写高阶微分方程

尸(乂-...””)/)=0

做变量代换处理：令,=>,必.，儿=得到

/•、

X(0)、'歹(0)、

y人(⑺v_%(0)八0)

Y'=•2，丫0-

*••*

J；1/(办<K(0);1尸(0),

例18求解描述振荡器的经典VerderP