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）
升级会员 