2.txt

这是一些关于学习matlab 的资料

q=dblquad(fun,xlower,xupper,ymin,ymax,tol,method,p1,p2,…)  %将可选参数p1,p2,..等传递给函数fun(x,y,p1,p2,…)。若tol=[]，method=[]，则使用缺省精度和算法。
例2-44
计算单位圆域上的积分： 
先把二重积分转化为二次积分的形式： 
f = inline(’exp(-x.^2/2).*sin(x.^2+y)’,’x’,’y’);
xlower = inline(’-sqrt(1-y.^2)’,’y’); xupper = inline(’sqrt(1-y.^2)’,’y’);
Q = quad2dggen(fun,xlower,xupper,-1,1,1e-4)
计算结果为：
  Q = 
      0.5368603818
2.4  常微分方程数值解
函数  ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb
功能  常微分方程（ODE）组初值问题的数值解
参数说明：
solver为命令ode45、ode23,ode113,ode15s,ode23s,ode23t,ode23tb之一。
Odefun 为显式常微分方程y’=f(t,y)，或为包含一混合矩阵的方程M(t,y)*y’=f(t,y)。命令ode23只能求解常数混合矩阵的问题；命令ode23t与ode15s可以求解奇异矩阵的问题。
Tspan 积分区间（即求解区间）的向量tspan=[t0,tf]。要获得问题在其他指定时间点t0,t1,t2,…上的解，则令tspan=[t0,t1,t2,…,tf]（要求是单调的）。
Y0 包含初始条件的向量。
Options 用命令odeset设置的可选积分参数。
P1,p2,… 传递给函数odefun的可选参数。
格式  [T,Y] = solver(odefun,tspan,y0)   %在区间tspan=[t0,tf]上，从t0到tf，用初始条件y0求解显式微分方程y’=f(t,y)。对于标量t与列向量y，函数f=odefun(t,y)必须返回一f(t,y)的列向量f。解矩阵Y中的每一行对应于返回的时间列向量T中的一个时间点。要获得问题在其他指定时间点t0,t1,t2,…上的解，则令tspan=[t0,t1,t2,…,tf]（要求是单调的）。
[T,Y] = solver(odefun,tspan,y0,options)   %用参数options（用命令odeset生成）设置的属性（代替了缺省的积分参数），再进行操作。常用的属性包括相对误差值RelTol（缺省值为1e-3）与绝对误差向量AbsTol（缺省值为每一元素为1e-6）。
[T,Y] =solver(odefun,tspan,y0,options,p1,p2…) 将参数p1,p2,p3,..等传递给函数odefun，再进行计算。若没有参数设置，则令options=[]。
1．求解具体ODE的基本过程：
（1）根据问题所属学科中的规律、定律、公式，用微分方程与初始条件进行描述。
F(y,y’,y’’,…,y(n),t) = 0
y(0)=y0,y’(0)=y1,…,y(n-1)(0)=yn-1
而y=[y;y(1);y(2);…,y(m-1)]，n与m可以不等
（2）运用数学中的变量替换：yn=y(n-1),yn-1=y(n-2),…,y2=y1=y，把高阶（大于2阶）的方程（组）写成一阶微分方程组： ， 
（3）根据（1）与（2）的结果，编写能计算导数的M-函数文件odefile。
（4）将文件odefile与初始条件传递给求解器Solver中的一个，运行后就可得到ODE的、在指定时间区间上的解列向量y（其中包含y及不同阶的导数）。
2．求解器Solver与方程组的关系表见表2-3。
表2-3
函数指令	含  义	函  数	含    义
求解器
Solver	ode23	普通2-3阶法解ODE	odefile	包含ODE的文件
	ode23s	低阶法解刚性ODE	选项	odeset	创建、更改Solver选项
	ode23t	解适度刚性ODE		odeget	读取Solver的设置值
	ode23tb	低阶法解刚性ODE	输出	odeplot	ODE的时间序列图
	ode45	普通4-5阶法解ODE		odephas2	ODE的二维相平面图
	ode15s	变阶法解刚性ODE		odephas3	ODE的三维相平面图
	ode113	普通变阶法解ODE		odeprint	在命令窗口输出结果
3．因为没有一种算法可以有效地解决所有的ODE问题，为此，MATLAB提供了多种求解器Solver，对于不同的ODE问题，采用不同的Solver。
表2-4  不同求解器Solver的特点
求解器Solver	ODE类型	特点	说明
ode45	非刚性	一步算法；4，5阶Runge-Kutta方程；累计截断误差达(△x)3	大部分场合的首选算法
ode23	非刚性	一步算法；2，3阶Runge-Kutta方程；累计截断误差达(△x)3	使用于精度较低的情形
ode113	非刚性	多步法；Adams算法；高低精度均可到10-3～10-6	计算时间比ode45短
ode23t	适度刚性	采用梯形算法	适度刚性情形
ode15s	刚性	多步法；Gear’s反向数值微分；精度中等	若ode45失效时，可尝试使用
ode23s	刚性	一步法；2阶Rosebrock算法；低精度	当精度较低时，计算时间比ode15s短
ode23tb	刚性	梯形算法；低精度	当精度较低时，计算时间比ode15s短
4．在计算过程中，用户可以对求解指令solver中的具体执行参数进行设置（如绝对误差、相对误差、步长等）。
表2-5  Solver中options的属性
属性名	取值	含义
AbsTol	有效值：正实数或向量
缺省值：1e-6	绝对误差对应于解向量中的所有元素；向量则分别对应于解向量中的每一分量
RelTol	有效值：正实数
缺省值：1e-3	相对误差对应于解向量中的所有元素。在每步(第k步)计算过程中，误差估计为：
e(k)<=max(RelTol*abs(y(k)),AbsTol(k))
NormControl	有效值：on、off
缺省值：off	为‘on’时，控制解向量范数的相对误差，使每步计算中，满足：norm(e)<=max(RelTol*norm(y),AbsTol)
Events	有效值：on、off	为‘on’时，返回相应的事件记录
OutputFcn	有效值：odeplot、odephas2、odephas3、odeprint
缺省值：odeplot	若无输出参量，则solver将执行下面操作之一：
画出解向量中各元素随时间的变化；
画出解向量中前两个分量构成的相平面图；
画出解向量中前三个分量构成的三维相空间图；
随计算过程，显示解向量
OutputSel	有效值：正整数向量
缺省值：[]	若不使用缺省设置，则OutputFcn所表现的是那些正整数指定的解向量中的分量的曲线或数据。若为缺省值时，则缺省地按上面情形进行操作
Refine	有效值：正整数k>1
缺省值：k = 1	若k>1，则增加每个积分步中的数据点记录，使解曲线更加的光滑
Jacobian	有效值：on、off
缺省值：off	若为‘on’时，返回相应的ode函数的Jacobi矩阵
Jpattern	有效值：on、off
缺省值：off	为‘on’时，返回相应的ode函数的稀疏Jacobi矩阵
Mass	有效值：none、M、
M(t)、M(t,y)
缺省值：none	M：不随时间变化的常数矩阵
M(t)：随时间变化的矩阵
M(t,y)：随时间、地点变化的矩阵
MaxStep	有效值：正实数
缺省值：tspans/10	最大积分步长
例2-45  求解描述振荡器的经典的Ver der Pol微分方程 
y(0)=1，y’(0)=0
令x1=y，x2=dy/dx，则
dx1/dt = x2
dx2/dt = μ(1-x2)-x1
编写函数文件verderpol.m：
function xprime = verderpol(t,x)
global MU
xprime = [x(2);MU*(1-x(1)^2)*x(2)-x(1)];
再在命令窗口中执行：
>>global MU
>>MU = 7;
>>Y0=[1;0]
>>[t,x] = ode45(‘verderpol’,0,40,Y0);
>>x1=x(:,1);x2=x(:,2);
>>plot(t,x1,t,x2)
图形结果为图2-20。
 
图2-20  Ver der Pol微分方程图
2.5  偏微分方程的数值解
MATLAB提供了一个专门用于求解偏微分方程的工具箱—PDE Toolbox (Paticial Difference Equation )。在本章，我们仅提供一些最简单、经典的偏微分方程，如：椭圆型、双曲型、抛物型等少数的偏微分方程，并给出求解方法。用户可以从中了解其解题基本方法，从而解决相类似的问题。
Matlab能解决的偏微分类型
    其中u = u(x,y)， 
 ， 
c = c(x,y) ， ， 
2.5.1  单的Poission方程
Poission方程是特殊的椭圆型方程： ， 
即 c = 1，a = 0，f = -1
Poission的解析解为： 。在下面计算中，用求得的数值解与精确解进行比较，看误差如何。
方程求解
（1）问题输入：
  c = 1；a = 0；f = 1；%方程的输入。给c，a，f赋值即可
  g = 'circleg'       %  区域G，内部已经定义为 circleg
   b = 'circleb1'      %  u在区域G的边界上的条件，内部已经定义好
（2）对单位圆进行网格化，对求解区域G作剖分，且作的是三角分划：
  [p,e,t] = initmesh（g,'hmax'，1）
（3）迭代求解：
  error = []；err = 1；
  while err > 0.001，
  [p,e,t]=refinemesh（'circleg',p,e,t）；
  u=assempde（'circleb1',p,e,t,1,0,1）；
  exact=-（p（1，：）^2+p（2，：）^2-1）/4；
  err=norm（u-exact',inf）；
  error=[error,err]；
  end
（4）误差显示与区域G内的剖分点数：
  Error: 1.292265e-002. Number of nodes: 25
  Error: 4.079923e-003. Number of nodes: 81
  Error: 1.221020e-003. Number of nodes: 289
  Error: 3.547924e-004. Number of nodes: 1089
（5）结果显示：
  subplot（2,2,1），pdemesh（p,e,t）%结果显示
  title（'数值解'）
  subplot（2,2,2），pdesurf（p,t,u）%精确解显示
  title（'精确解'）
  subplot（2,2,3），pdesurf（p,t,u-exact'）%与精确解的误差
  title（'计算误差'）
 
图2-21  Poission方程图
2.5.2  双曲型偏微分方程
1．Matlab能求解的类型： 
其中 ， ， ， ， ， 
2．形传递问题： ， 
即：c =1; a = 0; f = 0; d = 1
3．方程求解
（1）问题的输入：
  c = 1; a = 0; f = 0; d = 1;   % 输入方程的系数
  g = 'squareg'    % 输入方形区域G,内部已经定义好
  b = 'sqareb3'    % 输入边界条件，即初始条件
（2）对单位矩形G进行网格化：
  [p,e,t] = initmesh('squareg')
（3）定解条件和求解时间点：
  x = p(1,:)'; y = p(2,:)';
  u0 = atan(cos(pi/2*x));
  ut0 = 3*sin(pi*x).*exp(sin(pi/2*y));
  n = 31;
  tlist = linspace(0,5,n);
（4）求解：
  uu = hyperbolic(uo, ut0,tlist,b,p,e,t,c,a,f,d);
结果显示：计算过程中的时间点和信息
Time：0.166667
Time：0.333333
Time：4.33333
……
Time：4.66667
Time：4.83333
Time：5
428 successful steps
62 failed attempts
982 function evaluations
1 partial derivatives
142 LU decompositions
981 solutions of linear systems
（5）动画显示：
delta=-1:0.1:1;
[uxy,tn,a2,a3]=tri2grid(p,t,uu(:,1),delta,delta);
gp=[tn;a2;a3];
umax=max(max(uu));
umin=min(min(uu));
newplot;M=moviein(n);
for i=1:n,
    pdeplot(p,e,t,'xydata',uu(:,i),'zdata',uu(:,i),…
    'mesh','off','xygrid','on','gridparam',gp,…
    'colorbar','off','zstyle','continuous');
    axis([-1 1 -1 1 umin umax]);
    caxis([umin umax]);
    M(:,i)=getframe;
end
movie(M,5)
图2-22是动画过程中的一个状态。
 
图2-22  波动方程动画中的一个状态
2.5.3  抛物型偏微分方程
1．Matlab能求解的类型： 
其中 ， ， ， ， ， 
2．热传导方程： ， 
即：c = 1; a = 0; f = 1; d = 1;
3．问题计算
（1）问题的输入：
c = 1; a = 0; f = 1; d = 1;  %  输入方程的系数
g = 'squareg';  %  输入方形区域G
b = 'squareb1';  %  输入边界条件
（2）对单位矩形的网格化：
[p,e,t] = initmesh(g);
（3）定解条件和求解的时间点：
u0 = zeros(size(p, 2), 1);
ix = find(sqrt(p(1, :).^2+p(2, :).^2) < 0.4);
u0(ix) = ones(size(ix));
nframes = 20;
tlist=linspace(0,0.1,nframes)  % 在时间[0, 0.1]内20个点上计算，生成20帧
（4）求解方程：
u1 = parabolic(u0, tlist, b, p, e, t, c, a, f, d)
计算结果如下：
Time: 0.00526316
Time: 0.0105263
……
Time: 0.0947368
Time: 0.1
75 successful steps
1 failed attempts
154 function evaluations
1 partial derivatives
17 LU decompositions
153 solutions of linear systems
（5）动画显示：
x = linspace(-1,1,31); y = x;
newplot;
Mv = moviein(nframes);
umax=max(max(u1));
umin=min(min(u1));
for j=1:nframes
   u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));…
   surf(x,y,u);caxis([umin umax]);colormap(cool),…
   axis([-1 1 -1 1 0 1]);…
   Mv(:,j) = getframe;…
end
movie(Mv,10)
图2-23是动画过程中的瞬间状态。
 
图2-23  热传导方程动画瞬间状态图