遗迹之什么什么华秋实:MATLAB学习（序列1）偏微分方程组的求解 - blare-lover的日志 - 网易博...

MATLAB学习（序列1）偏微分方程组的求解

默认分类 2010-10-23 16:14:25 阅读51 评论0   字号：大中小 订阅

ode23    解非刚性微分方程，低精度，使用Runge-Kutta法的二三阶算法。
ode45    解非刚性微分方程，中等精度，使用Runge-Kutta法的四五阶算法。
ode113   解非刚性微分方程，变精度变阶次Adams-Bashforth-Moulton PECE算法。
ode23t    解中等刚性微分方程，使用自由内插法的梯形法则。
ode15s    解刚性微分方程，使用可变阶次的数值微分（NDFs）算法。
ode23s    解刚性微分方程，低阶方法，使用修正的Rosenbrock公式。
ode23tb    解刚性微分方程，低阶方法，使用TR-BDF2方法，即Runger-Kutta公式的第一级采用梯形法则，第二级采用Gear法。
[t，YY]=solver('F',tspan,Yo)
解算ODE初值问题的最简调用格式。
solver指上面的指令。

tspan=[0,30];     ％时域t的范围
y0=[1;0];         ％y（1）y(2)的初始值
[tt,yy]=ode45(@DyDt,tspan,y0);
plot(tt,yy(:,1)),title('x(t)')

function ydot=DyDt(t,y)
ydot=[y(2); 2*(1-y(1)^2)*y(2)-y(1)]

刚性方程：刚性是指其Jacobian矩阵的特征值相差十分悬殊。在解的性态上表现为，其中一些解变化缓慢，另一些变化快，且相差较悬殊，这类方程常常称为刚性方程，又称为Stiff方程。
刚性方程和非刚性方程对解法中步长选择的要求不同。
刚性方程一般不适合由ode45这类函数求解，而应该采用ode15s等。
如果不能分辨是否是刚性方程，先试用ode45，再用ode15s。

[t，YY，Te，Ye,Ie] = solver('F',tspan,Yo,options,p1,p2,…)
解算ODE初值问题的最完整调用格式。
为了能够解出方程，要用指令odeset确定求解的条件和要求。在MATLAB中，求解方程组的指令都有默认的求解的条件和要求（由结构数组options表示），但可以用odeset修改或重新建立，也可以用odeget去获取已有的“优化选项”的信息。指令odeset和odeget用法介绍如下：
语句格式如下：
options=odeset(‘name1’,value1,’name2’,value2,…)
options=odeset(oldopts,‘name1’,value1,’name2’,value2,…)
options=odeset(oldopts,newopts)
odeset
第一种调用格式是指定各个参数的取值，对不指定取值的参数，取默认值。在不引起混淆的情况下，参数名可以只键入前面的几个字母，也不必区分大小写，如用“abst”表示AbsTol.但数值的输入必须格式正确，否则仍采用默认值。
第二种格式使用了原来的优化选项，但对其中的参数1等指定了新值。
第三种格式合并了两个优化选项oldopts newopts，重复部分取newopts的指定值）：
第四种格式可在屏幕上显示如下全部可设置的参数及其默认值。

键入help odeset可查看全部参数的说明，下面对其中几个参数举例说明。
RelTol               相对误差，默认值为1e-3
AbsTol              绝对误差，默认值为1e-6
OutputFcn    输出方式，
默认值为‘odeplot’,其它选项有：
          odeplot       按时间顺序画出全部变量的解
          odephase2    二维相空间中两个变量的图形
          odephase2　三维相空间中三个变量的图形
          odeprint       打印输出
语句格式：
val=odeget(options,’name’)
这里options是由odeset设定的优化选项。
该语句从优化选项中提取指定的参数的取值。如果该参数没有指定值，则返回空阵［］。
options
odeget(options,'Reltol')
options=odeset(options,'Reltol',1e-6)
odeget(options,'Reltol')

function ydot=lorenfcn(t,y)
ydot=[-8/3*y(1)+y(2)*y(3);
    -10*y(2)+10*y(3);
    -y(2)*y(1)+35*y(2)-y(3)];

axis([10 50 -50 50 -50 50])
view(3)
hold on
title('Lorenz Attractor')
options=odeset('OutputFcn','odephas3');
[t,y]=ode23(@lorenfcn,[0 20],[0,0,eps],options);

％ EPS    Floating point relative accuracy.
    EPS returns the distance from 1.0 to the next largest
    floating point number. EPS is used as a default tolerance by
    PINV and RANK, as well as several other MATLAB functions.