matlab数值分析(2)

作者：未知  信息来源:未知  2006-1-26

字体大小:小 中 大网友评论条  进入论坛  

令y1=x 且 y2=dy/dx则 dy1/dt=y2根据这个微分方程组，可用MATLAB的函数ode23和ode45求出系统随时间变化的运动情况。调用这些函数时，需要编写一个函数M文件，给定当前时间及y1和y2的当前值，该函数返回上述导数值。所得函数M文件是：function yprime=vdpol(t , y)。 % output must be a column给定这个完整地描述微分方程的...

13.6  微分方程

一般微分方程式描述系统内部变量的变化率如何受系统内部变量和外部激励，如输入，的影响。当常微分方程式能够解析求解时，可用MATLAB的符号工具箱中的功能找到精确解。在本书的后面将介绍该工具箱的一些特点。

在微分方程难以获得解析解的情况下，可以方便地在数值上求解。为了说明起见，考虑描述振荡器的经典的范得波（Var der Pol）微分方程。

与所有的数值求解微分方程组的方法一样，高阶微分方程式必须等价地变换成一阶微分方程组。对于上述微分方程，通过重新定义两个新的变量，来实现种变换。

令y₁=x 且 y₂=dy/dx

则 dy1/dt=y₂

根据这个微分方程组，可用MATLAB的函数ode23和ode45求出系统随时间变化的运动情况。调用这些函数时，需要编写一个函数M文件，给定当前时间及y1和y2的当前值，该函数返回上述导数值。MATLAB中，这些导数由一个列向量给出。在本例中，这个列向量为yprime。同样，y1和y2合并写成列向量y。所得函数M文件：

function yprime=vdpol(t , y);

%  VDPOL(t , y) returns derivatives of the Van der Pol equation:

% 

%  x ‘‘-mu *(1-x ^2)*x ‘+x=0 (‘ = d/dx , ‘‘ = d^2/dx^2)

% 

%  let y(1)=x and y(2)=x‘

% 

%  then y(1) ‘ = y(2)

%      y(2) ‘ = MU*(1-y(1)^2)*y(2)-y(1)

global MU %   choose 0<MU<10 in Command workspace

yprime=[y(2)  MU*(1-y(1)^2)*y(2)-y(1)];  %  output must be a column

给定这个完整地描述微分方程的函数，计算结果如下：

>>global MU %  define MU as a global variable in the Command Workspace

>>MU=2;  %  set global parameter to desired value

>>[t , y]=ode23(‘ vdpol ‘ , 0 , 30 , [1 ; 0]);  %  to=0 , tf=30 , yo=[1 ; 0]

>>y1=y( : , 1);  %  first column is y(1) versus time points in t

>>y2=y( : , 2);  %  second column is y(2)

>>plot(t , y1 , t , y2 , ‘ -- ‘)

>>xlabel(‘ Time , Second ‘) , ylabel(‘ Y(1) and Y(2) ‘)

>>title(‘ Van der Pol Solution for mu=2 ‘)

所得的图见图13.9。

图13.9当mu=2时的范得波方程的运动曲线

在图13.9中，y2（虚线）是y1（实线）的导数。传递给ode23的数由ode23(f_name , to , tf , yo , to1)描述。这里f_name是计算导数的M文件函数的字符串名，to是初始时间，tf是终止时间，yo是初始条件向量。可选择的参数to1（缺省值to1=1e-3）是所需的相对精度。在上例中，起始时间是第0秒，终止时间是第30秒，初始条件为y=[1;0]。两个输出参数是列向量t和矩阵y，其中向量t包含了估计响应的时间点，而矩阵y的列数等于微分方程组的个数（本例为2），且其行数与t相同。t中的时间点不是等间隔的，因为为了保持所需的相对精度，积分算法改变了步。

函数ode45的使用与ode23完全一样。两个函数的差别在于必须与所用的内部算法相关。两个函数都运用了基本的龙格-库塔(Runge-Kutta)数值积分法的变形。ode23运用一个组合得2/3阶龙格-库塔-芬尔格(Runge-Kutta-Fehlerg)算法，而ode45运用组合的4/5阶龙格-库塔-芬尔格算法。一般地，ode45可取较多的时间步。所以，要保持与ode23相同误差时，在to和tf之间可取较少的时间步。然而，在同一时间，ode23每时间步至少调用f_name 3次，而ode45每时间步至少调用f_name 6次。

正如使用高阶多项式内插常常得不到最好的结果一样，ode45也不总是比ode23好。如果ode45产生的结果，对作图间隔太大，则必须在更细的时间区间，对数据进行内插，比如用函数interp1。这个附加时间点会使ode23更有效。作为一条普遍规则，在所计算的导数中，如有重复的不连续点，为保持精度致使高阶算法减少时间步长，这时低阶算法更有效。正是由于这个原因，电子电路分析按缺省，就用一阶算法编程，并且最多提供二阶算法来解决暂态时间响应问题。此外，通过对tol设置更小的值，要达到更高的精度，没有必要使绝对误差更小。tol设置每时间步的相对精度，不一定引起绝对误差减少。

总之，不要盲目使用数值方法。对于给定的问题，在决定最好的方法之前，要试一试各种可能的方法。有关微分方程数值解法的更进一步信息，请参考数值分析方面的书籍。有些参考书还提供了一些关于算法选择和如何处理那些时间常数变化范围大的病态方程的非

 

13.7  M文件举例

这里所介绍的《精通MATLAB工具箱》中的M文件可近似求解由采样值给出的函数的积分和微分。这里假定这些函数本身不存在，且独立变量也许不是线性间隔。例如，已装载到MATLAB中要分析的数据来源于实验测试。

对于所包含的数据缺乏函数描述，有许多种积分和微分的方法。如前所述，们可以用最小二乘多项式拟合数据，然后在多项式的描述上进行操作。另一种方法是寻找数据的三次样条表示，然后运用《精通MATLAB工具箱》中的函数spintgrl和spderiv来分别寻找积分和微分的样条表示。这里所介绍的方法提供了另一种更简单的方法。积分用梯形规则计算。用加权中心差分计算微分。此外，将函数设计成在矩阵形式下工作，矩阵的列代表各与自变量有关的因变量。

正如这章前面所述，MATLAB函数trapz计算在某有限区间的梯形积分。这里我们寻找的积分是自变量为x的函数。即如果y=f(x)，我们寻找：

式中的x1是向量x的第一个元素。用梯形规则，这个积分近似为：

  且S(x₁)=0

这样，第k个数据点的积分是上述梯形面积的累加和。函数mmintgrl实现的这个算法如下：

function z=mmintgrl(x , y)

%  MMINTgrl Compute Integral using Trapezoidal Rule.

%   MMINTGRL(X , Y) computes the integral of the function y=f(x) given the

%   data in X and Y. X must be a vector ,  but Y may be a column oriented

%   data matrix. The length of X must equal the length of Y if Y is a

%   vector ,  or it must equal the number of rows in Y if Y is a matrix.

% 

%   X need not be equally spaced. The trapezoidal algorithm is used.

% 

%   See also mmderiv

%   Copyrigth (c) 1996 by Prentice-Hall , Inc.

flag=0;                        %   falg is True if y is a row

x=x( : ); nx=length(x);  %   make x a column

[ry , cy]=size(y);

if ry==1&cy==nx , y=y .' ; ry=cy ; cy=1 ; flag=1 ; end

if nx~=ry ,  error(' X and Y not the right size ') , end

dx=x(2 : nx)-x(1 : nx-1);             %   width of each trapezoid

dx=dx( : , ones(1 , cy));              %   duplicate for each column in y

yave=(y(2 : ry , : )+y(1 : ry-1 , :))/2;                 %   average of heights

z=[zeros(1 , cy); cumsum(dx .* yave)];            %   Use cumsum to find area

if flag , z=z'; end                                      %   if y was a row , return a row

在介绍上述函数的使用之前，考虑微分。在这种情况下，人们感兴趣的就是刚给定数据点的近似斜率。这里介绍一种下述的中心差分

图13.10 加权中心差分方法

从图13.11可知，在第k个点的近似微分是：

       式中，

并且M_k是连接y_k-1到y_k的直线的斜率。这样，第k点的微分是相邻两点间斜率的加权平均，离该点越近的点权越重。在第一个和最后一个数据点上，不能简单按照上述方法进行处理，因为这两个点都没有伴随的直线段。对于这些数据点，需要用另外的方法。这里所采取的方法是用二次多项式拟合前3个点（或最后3个点），并且计算这个多项式第一个（或最后一个）点的微分。函数mmderiv实现的这个算法如下：

function z=mmderiv(x , y)

%  MMDERIV Compute Derivative Using Weighted Central Differences.

%   MMDERIV(X , Y) computes the derivative of the function y=f(x) given the

%   data in X and Y. X must be a vector , but Y may be a column oriented

%   data matrix. The length of X must equal the length of Y if Y is a

%   vector ,  or it must equal the number of rows in Y if Y is a matrix.

% 

%   X need not be equally spaced.Weighted central difference are used.

%   Quadratic approximation is used at the endpoints.

%  

%   See also mmintgrl

%   Copyrigth (c) 1996 by Prentice-Hall , Inc.

flag=0;                        %   flag is True if y is a row

x=x( : ); nx=length(x);  %   make x a column

[ry , cy]=size(y);

if ry==1&cy==nx , y=y .'; ry=cy; cy=1; flag=1; end

if nx~=ry ,  error(' X and Y not the right size ') , end

if nx<3 ,  error(' X and Y must have st least three elements ') , end

dx=x(2 : nx)-x(1 : nx-1);                    %   first difference in x

dx=dx+(dx==0)*eps;                  %   make infinite slopes finite

dxx=x(3 : nx)-x(1 : nx-2);          %   second difference in x

dxx=dxx+(dxx==0)*eps;             %   make infinite slopes finite

alpha=dx(1 : nx-2) ./ dxx                    %   central difference weight

alpha=alpha( : , ones(1 , cy));             %   duplicate for each column in y

dy=y(2 : ry , :)-y(1 : ry-1 , : );                   %   first difference in y

dx=dx( : , ones(1 , cy));                     %   duplicate dx for each column in y

%   now apply weighting to dy

z=alpha .* dy(2:ry-1 , :) ./ dx(2 : nx-1 , : )+(1-alpha) .* dy(1 : ry-2 , : ) ./ dx(1 : nx-2 , : );

z1=zeros(1 , cy)>=z1;

for i=1 : cy    %   fit quadratic at endpoints of each column

   p1=polyfit(x(1 : 3) , y(1 : 3 , i) , 2);                    %   quadratic at first point

   z1(i)=2*p1(1)*x(1)+p1(2);                   %   evalute poly derivative

   pn=polyfit(x(nx-2 : nx) , y(ry-2 : ry , i) , 2);%   quadratic at last point   

   zn(i)=2*pn(1)*x(nx)+pn(2);                 %   evaluate poly derivative

end

z=[z1;  z;  zn];

if flag , z=z'; end                 %   if y was a row , return a row

最后，给出一个例子：

>>x=linspace(0 , 2*pi , 30)

>>y=sin(x);    %  create data

>>yi=mmintgrl(x , y);  %  find integral

>>yd=mmderiv(x , y);  %  find derivative

>>plot(x , y , x , yi , ‘ - ‘ , x , yd , ‘ : ‘)   %  plot results

注意这个积分定性地证明了等式：

而微分定性地证明了等式：

图13.11 y=sin(x)极其积分、微分曲线

13.8  小结

表13.1总结了本章所讨论的函数。

表13.1