Cl_matcont是一个常微分方程(组)的MatLab分岔软件包,其中有个system文件夹,在这里可以定义方程或
方程组,但其中有几个矩阵如jacobian,jacobianp,hessian,hessiansp,der3,der4,der5,一般资料中没有给出怎样计算它们的方法,如果算不出这些矩阵就没法进行分岔分析模拟,现将这几个矩阵的算法注解如下,即%后的注解。 假设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;
|