Question

高手帮忙~~~帮忙分析一段MATLAB程序~~~~高加分！！！

帮忙分析一段MATLAB的程序，大概说说每一步都是什么意思，挺着急的。。。。救急。。。。保证还有高加分！！！！！感激不尽。。。。高手帮忙。。。

function OC_118_569_1995_fig1
clear all;
close all;
clc;
clf;

r=[0.1 1 10 40];
beta=35.3;
options=odeset('RelTol',1e-6,'AbsTol',1e-6,'Stats','off','Events',@events);
xmax=2;
y0=[1;0];
n=1;
while n<=4
[Xr,Yr]=ode45(@func,[0 xmax],y0,options,beta,r(n));
[Xl,Yl]=ode45(@func,[0 -xmax],y0,options,beta,r(n));
X=[flipud(Xl(2:end));Xr];
Y=[flipud(Yl(2:end,1));Yr(:,1)];

if n==1
plot(X,Y.^2,'k:','Linewidth',1.5);hold on;
end
if n==2
plot(X,Y.^2,'k-.','Linewidth',1.5);hold on;
end
if n==3
plot(X,Y.^2,'k','Linewidth',1.5);hold on;
end
if n==4
plot(X,Y.^2,'k--','Linewidth',1.5);hold on;
end

n=n+1;
end
legend('r=0.1','r=1','r=10','r=40',2);%0:默认,1:右上侧,2:左上侧,3:左下侧,4:右下侧,-1:右侧外边
legend('boxoff');
set(gca,'xticklabel',[-2:.5:2],'FontSize',14);
set(gca,'yticklabel',[0:.1:1],'FontSize',14);
% set(gca,'xgrid','on');
xlabel('\it s','Fontsize',20);
ylabel('Normalized Intensity','Fontsize',20);
text(1,0.8,'{\it\beta} =35.3','Fontsize',20);

%==========================================================================
% 控制ode45的事件函数
function [value,isterminal,direction]=events(x,y,beta,r)
value=[y(1);y(2)];
isterminal=[1;1];
direction=[0;0];
%==========================================================================
% 微分方程(5)
function dy=func(x,y,beta,r)
dy=[y(2);
-2*beta*y(1)*log(1+r)/r+2*beta*y(1)/(1+r*y(1)^2)];

做出来的图我插入了，请高手帮忙读下程序。。。

Answer 1

常微分（ODE）方程的数值求解器有：非刚性求解器(计算的精度从低到高)ode23，ode45，ode113，刚性方程求解器(适用的刚性从弱到强) ode15s，ode23s，ode23t，ode23tb，隐方程求解器ode15i
所谓刚性方程，就是指它的解的曲线有剧烈的或缓慢的变化。如van der Pol方程(教材p144例4)就是一个刚性方程.
Ode求解器默认的相对误差是1e-3，绝对误差是1e-6，要改变默认的精度设置，可以用odeset来设定Options。具体设置方法，求助于help功能。Ode求解器中可以求解带有参数的微分方程。
常用的精度设置如
Options=odeset(‘RelTol’, 1e-5,’AbsTol’, [1e-8, 1e-7]);
其中绝对误差可以对每个未知函数的分量分别规定，写成一个向量，维数等于方程的维数。如对各分量的绝对误差设置相同则只须写一个标量误差。
每一积分步第i个分量的误差满足e(i) <= max(RelTol*abs(y(i)),AbsTol(i)).
如果只要对解的范数作误差控制，而不需对解的每个分量作误差控制，则在Options中可以加上选项’NormControl’, ‘on’. 这时每一积分步误差的范数满足norm(e) <= max(RelTol*norm(y),AbsTol). 这个选项对那些解的范数会等于零的方程特别有用，不用此选项时，为了要达到苛刻的误差要求，步长会取得很小，将大大减慢求解过程以致求解失败.
时滞是常数的时滞常微分方程DDE的数值求解器有dde23，要改变默认的设置，可以用ddeset来设定。

常微分方程求解器的options还可以设置一个有用的功能，语法是
Options=odeset( 'Events',@EVENTS);(Options可以是自定义名，其中各种设置如精度设置可以写在一起，用逗号分开)，用ode45求解时，可用格式
[TOUT,YOUT,TE,YE,IE] = ode45(@ODEFUN,TSPAN,Y0,Options, P)
在odeset中设置了一个事件函数@EVENTS，是自编m函数，函数名自定。函数的格式是
[VALUE,ISTERMINAL,DIRECTION] = EVENTS(T,Y,P).
事件函数的输入是和微分方程的函数输入相同，顺序相同. 事件函数的输出是3个列向量. 例如 [VALUE,ISTERMINAL,DIRECTION]，名称自定，列向量的维数是事件的个数，VALUE(I) 是事件函数的第I个分量表达式的值, ISTERMINAL(I)=1表示事件函数的第I个分量的值等于零时积分终止，不然等于0 ， DIRECTION(I)=0 表示要计算事件函数的第I个分量的所有的零点 (默认), +1 仅计算事件函数第I个分量在零点递增的零点 -1仅计算事件函数的第I个分量在零点递减的零点.
TSPAN是微分方程求解区间，例如可用 [0 pi]表示，也可以指定TOUT为求解区间中的一些点列，如TSPAN=0: 0.1: 3.2; Y0是微分方程初值问题中的初值，如Y0=[0 0]；
TOUT是微分方程求解区间TSPAN中的时间点列，YOUT是时间点列TOUT上方程解的值，是矩阵形式，列数等于方程的维数，TOUT(:,J)是解的第J个分量在时间点列TOUT上的向量值。 TE是列向量，是事件发生的时刻序列. YE的各行向量是事件发生各时刻方程的解的向量值, IE表示在TE时刻发生的事件在事件函数中的序号.