分岔软件Cl_matcont中关于jacobian，jacobianp，hessian，hessiansp，der3矩阵的注解 转载自sinaBlog

  Cl_matcont是一个常微分方程（组）的MatLab分岔软件包，其中有个system文件夹，在这里可以定义方程或
方程组，但其中有几个矩阵如jacobian，jacobianp，hessian，hessiansp，der3,der4,der5，一般资料中没有给出怎样计算它们的方法，如果算不出这些矩阵就没法进行分岔分析模拟，现将这几个矩阵的算法注解如下，即%后的注解。（Cl_matcont下载地址为http://www.matcont.ugent.be/）

    假设system中有个bratu文件，编程如下,则其中的矩阵，如jacobian，jacobianp，hessian，hessiansp，der3,der4,der5等矩阵可作如下注解。

function out = bratu
out{1} = @init;
out{2} = @fun_eval;
out{3} = @jacobian;
out{4} = @jacobianp;
out{5} = @hessians;
out{6} = @hessiansp;
out{7} = @der3;
out{8} = [];
out{9} = [];
out{10}= @userf1;

% --------------------------------------------------------------------------
function dydt = fun_eval_r(t,kmrgd,a) %定义微分方程组
dydt=[-2*kmrgd(1)+kmrgd(2)+a*exp(kmrgd(1));;
kmrgd(1)-2*kmrgd(2)+a*exp(kmrgd(2));;];

% --------------------------------------------------------------------------
function [tspan,y0,options] = init
y0=[0,0];  %初始条件
options = odeset('Jacobian',handles(3),'JacobianP',handles(4),'Hessians',handles(5),'HessiansP',handles(6));
handles = feval_r(bratu);
tspan = [0 10];  %运行时间

% --------------------------------------------------------------------------
function jac = jacobian(t,kmrgd,a)
jac=[ a*exp(kmrgd(1)) - 2 , 1 ; 1 , a*exp(kmrgd(2)) - 2 ]; %dydt关于kmrgd(1)和kmrgd(2)计算雅可比矩阵
% --------------------------------------------------------------------------
function jacp = jacobianp(t,kmrgd,a)
jacp=[ exp(kmrgd(1)) ; exp(kmrgd(2)) ]; %dydt关于参数a计算雅可比矩阵
% --------------------------------------------------------------------------
function hess = hessians(t,kmrgd,a)
hess1=[ a*exp(kmrgd(1)) , 0 ; 0 , 0 ];
hess2=[ 0 , 0 ; 0 , a*exp(kmrgd(2)) ];  %dydt关于kmrgd(1)和kmrgd(2)计算hessian矩阵
hess(:,:,1) =hess1;
hess(:,:,2) =hess2;
% --------------------------------------------------------------------------
function hessp = hessiansp(t,kmrgd,a)
hessp1=[ exp(kmrgd(1)) , 0 ; 0 , exp(kmrgd(2)) ]; %雅可比矩阵关于参数a计算导数
hessp(:,:,1) =hessp1;
%---------------------------------------------------------------------------
function tens3  = der3(t,kmrgd,a) %共有2*2个导数矩阵，即变量数为2，则为2*2个导数矩阵
tens31=[ a*exp(kmrgd(1)) , 0 ; 0 , 0 ]; %hessian矩阵hess1关于kmrgd(1)求导数
tens32=[ 0 , 0 ; 0 , 0 ];  %hessian矩阵hess1关于kmrgd(2)求导数
tens33=[ 0 , 0 ; 0 , 0 ];  %hessian矩阵hess2关于kmrgd(1)求导数
tens34=[ 0 , 0 ; 0 , a*exp(kmrgd(2)) ]; %hessian矩阵hess2关于kmrgd(2)求导数
tens3(:,:,1,1) =tens31;
tens3(:,:,1,2) =tens32;
tens3(:,:,2,1) =tens33;
tens3(:,:,2,2) =tens34;
%---------------------------------------------------------------------------
function tens4  = der4(t,kmrgd,a) %对der3关于各变量再次求导，共2*2*2个矩阵
%---------------------------------------------------------------------------
function tens5  = der5(t,kmrgd,a) %对der4关于各变量再次求导，共2*2*2*2个矩阵

function userfun1=userf1(t,kmrgd,a)
 userfun1=a-0.2;