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

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

《Matlab笔记数值计算高数篇015.docx》由会员分享，可在线阅读，更多相关《Matlab笔记数值计算高数篇015.docx（13页珍藏版）》请在冰豆网上搜索。

Matlab笔记数值计算高数篇015

15.数值计算—高数篇

一、求极限

limit（f,x,a）——求极限

limit（f,x,a,'right'）——求右极限

limit（f,x,a,'left'）——求左极限

例1求

代码：

symsx;

y=（3*x^2+5）/（5*x+3）*sin（2/x）;

limit（y,x,inf）

运行结果：

ans=6/5

注：

Matlab求二元函数的极限，是用嵌套limit函数实现的，相当于求的是累次极限，需要注意：

有时候累次极限并不等于极限。

例2求

代码：

symsxab;

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阶导数，并绘图

代码：

symsxab;

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设

，求

代码：

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）^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.00001.0000

f1=2.7969e-010

四、求不定积分与定积分

1.符号积分

int（f,x）——求f（x）关于x的不定积分

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

例7求积分

和

代码：

symsxa;

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求二重积分

先化为累次积分：

原式=

代码：

symsxy;

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阶泰勒式，

在

的傅里叶展开。

代码：

symsax;

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求级数

symsn;

symsum（（-1）^（n+1）/（n+1）^2,n,1,inf）

运行结果：

ans=1-pi^2/12

七、代数方程

1.求代数方程的解析解

solve（‘eq1’,’eq2’,…,’var1’,’var2’,…）

例13解方程

的x和b，以及方程组

代码：

symsabcx;

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/3y=-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.

求解器

ODE类型

特点

说明

ode45

非刚性

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

方法累积截断误差

大部分场合的首选算法

ode23

非刚性

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

方法累积截断误差

使用于精度较低的情形

ode113

非刚性

多步法，Adams算法，高低精度均可达到

计算时间比ode45短

ode23t

适度刚性

采用梯形算法

适度刚性情形

ode15s

刚性

多步法，Gear’s反向

数值积分，精度中等

若ode45失效时，

可尝试使用

ode23s

刚性

一步法，2阶Rosebrock算法，

低精度。

当精度较低时，

计算时间比ode15s短

调用格式：

[T,X]=solver（odefun,tspan,X0）

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

做变量代换处理：

令

，得到

，

例18求解描述振荡器的经典VerderPol微分方程（取

）：

，

，

做变量代换处理，

，则

代码：

先编写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）