Initialize

> s:=rand(0..10000);

> Digits:=5;

> m:=3;

> n:=2;

> g:=matrix(n,n,[][]);

> g_inv:=matrix(n,n,[][]);

> alpha:=0.5;

> sub_alpha:=2/(1-alpha);

> add_alpha:=(1+alpha)/2;

> F:=matrix(2,3,[[0.8,0.5,0.1],[0.1,0.5,0.8]]);

Set Functions

> m_alpha:=(x,theta1,theta2)->((F[1,x]*theta1+F[2,x]*theta2)/sub_alpha)^sub_alpha;

> eta:=(j,theta1,theta2)->sum(F[j,'x']*m_alpha('x',theta1,theta2)^add_alpha, 'x'=1..m)/add_alpha;

> eta_q:=j->sum(F[j,'x']*(q['x']^add_alpha), 'x'=1..m)/add_alpha;

> psi:=(theta1,theta2)->sum(m_alpha('x',theta1,theta2),'x'=1..m)/add_alpha;

> d:=(tq1,tq2)->psi(tq1,tq2)-(theta1_q*eta_q(1)+theta2_q*eta_q(2));

> D_alpha:=(theta1,theta2)->psi(theta1,theta2)-(theta1*eta_q(1)+theta2*eta_q(2))-d(theta1_q,theta2_q);

> delta:=(theta1,theta2,theta3)->theta3-D_alpha(theta1,theta2);

Calculation g

> for i from 1 to n do for j from 1 to n do g[i,j]:=sum(F[i,'x']*F[j,'x']*m_alpha('x',theta1,theta2)^alpha,'x'=1..3) end do end do;

Calculation g Inverse

> inv:=proc(t::algebraic)local i;

> i:=t[1,1]*t[2,2]-t[1,2]*t[2,1];

> matrix(2,2,[[t[2,2]/i,-t[1,2]/i],[-t[2,1]/i,t[1,1]/i]]);

> end proc;

> g_inv:=inv(g);

> DD:=(j,theta1,theta2)->eta(j,theta1,theta2)-eta_q(j);

> q := [s()*0.001, s()*0.001,s()*0.001];

> RD:=(j,theta1,theta2)->sum(F[j,'x']*(m_alpha('x',theta1,theta2)^add_alpha-q['x']^add_alpha), 'x'=1..m)/add_alpha;

> init_theta:=[s()*0.0001,s()*0.0001];

> ode1:=-(g_inv[1,1]*RD(1,theta1,theta2)+g_inv[1,2]*RD(2,theta1,theta2));

> ode2:=-(g_inv[2,1]*RD(1,theta1,theta2)+g_inv[2,2]*RD(2,theta1,theta2));

> ff:=dsolve({diff(theta1(t),t)=ode1,diff(theta2(t),t)=ode2,theta1(0)=init_theta[1],theta2(0)=init_theta[2]},{theta1(t),theta2(t)}, type=numeric, output=listprocedure);

> with(plots):

> odeplot(ff,[[t,theta1(t)],[t,theta2(t)]],0..30,view=[0..30,-30..30],title=cat("alpha=",convert(alpha,string),"\nq=",convert(q,string),"\nInit_Theta=",convert(init_theta,string)),legend=["Theta1", "Theta2"]);

> with(plots):

> fieldplot([ode1,ode2],theta1=2..30,theta2=2..30,title=cat("alpha=",convert(alpha,string),"\nq=",convert(q,string),"\nInit_Theta=",convert(init_theta,string)));

> with(plots):

> implicitplot({ode1=0,ode2=0},theta1=1..30,theta2=1..30,title=cat("alpha=",convert(alpha,string),"\nq=",convert(q,string),"\nInit_Theta=",convert(init_theta,string)));

> conv:=fsolve({ode1=0,ode2=0},{theta1=1..30,theta2=1..30});

> theta1_q:=eval(theta1,conv);

> theta2_q:=eval(theta2,conv);

> tmp:=(i,j)->RD(i,theta1,theta2)*g_inv[i,j]*RD(j,theta1,theta2);

> ode3:=-delta(theta1,theta2,theta3)*(delta(theta1,theta2,theta3)+tmp(1,1)+tmp(1,2)+tmp(2,1)+tmp(2,2));

> ode4:=-delta(theta1,theta2,theta3)*(g_inv[1,1]*RD(1,theta1,theta2)+g_inv[1,2]*RD(2,theta1,theta2));

> ode5:=-delta(theta1,theta2,theta3)*(g_inv[2,1]*RD(1,theta1,theta2)+g_inv[2,2]*RD(2,theta1,theta2));

> ff2:=dsolve({diff(theta1(t),t)=ode4,diff(theta2(t),t)=ode5,diff(theta3(t),t)=ode3,theta1(0)=init_theta[1],theta2(0)=init_theta[2],theta3(0)=26},{theta1(t),theta2(t),theta3(t)}, type=numeric, output=listprocedure);

> with(plots):

> odeplot(ff2,[[t,theta1(t)],[t,theta2(t)],[t,theta3(t)]],0..5,view=[0..100,-100..100],title=cat("alpha=",convert(alpha,string),"\nq=",convert(q,string),"\nInit_Theta=",convert(init_theta,string),"\nconv Theta=[",convert(theta1_q,string),",",convert(theta2_q,string),"]"),legend=["Theta1", "Theta2","Theta3"]);

>

>