Cdc14 modell

Előkészítés: 4-es modellben a felírt egyenletek és a referencia-paraméterértékek:

Needs["PlotLegends`"]

valtozok={x[t],y[t],z[t],w[t]};

modeleq={x'[t]==a1-(a21+a22*y[t])*x[t],
y'[t]==(((a31+a32*(t4-w[t]))*(t1-y[t]))/(j3+(t1-y[t])))-(((a41+a42*x[t])*y[t])/(j3+y[t])),
z'[t]==(((a51*(t4-w[t])+a52*t2)*(t3-z[t]-w[t]))/(j6+(t3-z[t]-w[t])))-((a7*x[t]*z[t])/(j6+z[t]))-l1*z[t]*(t4-w[t])+l2*w[t],
w'[t]==l1*z[t]*(t4-w[t])-l2*w[t]-((a7*x[t]*w[t])/(j6+w[t]))};

inipar={x[0]==0.04,y[0]==0.33,z[0]==0.499,w[0]==0.490};

params={a1,a21,a31,a41,a51,a7,a22,a32,a42,a52,l1,l2,j3,j6,t1,t2,t3,t4} ;

refparam={a1->0.03,a21->0.03,a31->0.01,a41->0.05,a51->0.001,a7->2,a22->2,a32->1,a42->0.5,a52->0.45,l1->1000,l2->10,j3->0.001,j6->0.02,t1->1,t2->1,t3->1,t4->0.5};

nds=NDSolve[Join[modeleq,inipar]/.refparam,{x[t],y[t],z[t],w[t]},{t,0,800},MaxSteps->100000];

plot1=Plot[{x[t]/.nds,y[t]/.nds,z[t]/.nds,w[t]/.nds},{t,0,60},PlotStyle->{Red,Blue,Black,Green}, PlotLegend->{"x(t)","y(t)","z(t)","w(t)"},LegendPosition->{1,-0.5},ImageSize->400]

FindMinimum[x[t]/.nds,{t,14}]

FindMinimum[x[t]/.nds,{t,32}]

tau=32.8202-16.0856

A nyers érzékenység:

J[t]=D[#,{valtozok}]&/@modeleq[[All,2]];

lista1=First/@refparam

lista2=Transpose[{{xa1[t],ya1[t],za1[t],wa1[t]},{xa21[t],ya21[t],za21[t],wa21[t]},{xa31[t],ya31[t],za31[t],wa31[t]},{xa41[t],ya41[t],za41[t],wa41[t]},{xa51[t],ya51[t],za51[t],wa51[t]},{xa7[t],ya7[t],za7[t],wa7[t]},{xa22[t],ya22[t],za22[t],wa22[t]},{xa32[t],ya32[t],za32[t],wa32[t]},{xa42[t],ya42[t],za42[t],wa42[t]},{xa52[t],ya52[t],za52[t],wa52[t]},{xl1[t],yl1[t],zl1[t],wl1[t]},{xl2[t],yl2[t],zl2[t],wl2[t]},{xj3[t],yj3[t],zj3[t],wj3[t]},{xj6[t],yj6[t],zj6[t],wj6[t]},{xt1[t],yt1[t],zt1[t],wt1[t]},{xt2[t],yt2[t],zt2[t],wt2[t]},{xt3[t],yt3[t],zt3[t],wt3[t]},{xt4[t],yt4[t],zt4[t],wt4[t]}}];

Sjobbo[t]=J[t].lista2+Transpose[Map[{D[modeleq[[All,2]][[1]],lista1[[#]]],D[modeleq[[All,2]][[2]],lista1[[#]]],D[modeleq[[All,2]][[3]],lista1[[#]]],D[modeleq[[All,2]][[4]],lista1[[#]]]}&,Range[18]]];

Sbalo[t]=Transpose[{{xa1'[t],ya1'[t],za1'[t],wa1'[t]},{xa21'[t],ya21'[t],za21'[t],wa21'[t]},{xa31'[t],ya31'[t],za31'[t],wa31'[t]},{xa41'[t],ya41'[t],za41'[t],wa41'[t]},{xa51'[t],ya51'[t],za51'[t],wa51'[t]},{xa7'[t],ya7'[t],za7'[t],wa7'[t]},{xa22'[t],ya22'[t],za22'[t],wa22'[t]},{xa32'[t],ya32'[t],za32'[t],wa32'[t]},{xa42'[t],ya42'[t],za42'[t],wa42'[t]},{xa52'[t],ya52'[t],za52'[t],wa52'[t]},{xl1'[t],yl1'[t],zl1'[t],wl1'[t]},{xl2'[t],yl2'[t],zl2'[t],wl2'[t]},{xj3'[t],yj3'[t],zj3'[t],wj3'[t]},{xj6'[t],yj6'[t],zj6'[t],wj6'[t]},{xt1'[t],yt1'[t],zt1'[t],wt1'[t]},{xt2'[t],yt2'[t],zt2'[t],wt2'[t]},{xt3'[t],yt3'[t],zt3'[t],wt3'[t]},{xt4'[t],yt4'[t],zt4'[t],wt4'[t]}}];

erzekeny=Map[Thread[Sjobbo[t][[#]]==Sbalo[t][[#]]]&,Range[4]];

inierzekeny={xa1[0]==0,ya1[0]==0,za1[0]==0,wa1[0]==0,xa21[0]==0,ya21[0]==0,za21[0]==0,wa21[0]==0,xa31[0]==0,ya31[0]==0,za31[0]==0,wa31[0]==0,xa41[0]==0,ya41[0]==0,za41[0]==0,wa41[0]==0,xa51[0]==0,ya51[0]==0,za51[0]==0,wa51[0]==0,xa7[0]==0,ya7[0]==0,za7[0]==0,wa7[0]==0,xa22[0]==0,ya22[0]==0,za22[0]==0,wa22[0]==0,xa32[0]==0,ya32[0]==0,za32[0]==0,wa32[0]==0,xa42[0]==0,ya42[0]==0,za42[0]==0,wa42[0]==0,xa52[0]==0,ya52[0]==0,za52[0]==0,wa52[0]==0,xl1[0]==0,yl1[0]==0,zl1[0]==0,wl1[0]==0,xl2[0]==0,yl2[0]==0,zl2[0]==0,wl2[0]==0,xj3[0]==0,yj3[0]==0,zj3[0]==0,wj3[0]==0,xj6[0]==0,yj6[0]==0,zj6[0]==0,wj6[0]==0,xt1[0]==0,yt1[0]==0,zt1[0]==0,wt1[0]==0,xt2[0]==0,yt2[0]==0,zt2[0]==0,wt2[0]==0,xt3[0]==0,yt3[0]==0,zt3[0]==0,wt3[0]==0,xt4[0]==0,yt4[0]==0,zt4[0]==0,wt4[0]==0};

mo=NDSolve[Join[modeleq,inipar,erzekeny,inierzekeny]/.refparam,{x[t],y[t],z[t],w[t],xa1[t],ya1[t],za1[t],wa1[t],xa21[t],ya21[t],za21[t],wa21[t],xa31[t],ya31[t],za31[t],wa31[t],xa41[t],ya41[t],za41[t],wa41[t],xa51[t],ya51[t],za51[t],wa51[t],xa7[t],ya7[t],za7[t],wa7[t],xa22[t],ya22[t],za22[t],wa22[t],xa32[t],ya32[t],za32[t],wa32[t],xa42[t],ya42[t],za42[t],wa42[t],xa52[t],ya52[t],za52[t],wa52[t],xl1[t],yl1[t],zl1[t],wl1[t],xl2[t],yl2[t],zl2[t],wl2[t],xj3[t],yj3[t],zj3[t],wj3[t],xj6[t],yj6[t],zj6[t],wj6[t],xt1[t],yt1[t],zt1[t],wt1[t],xt2[t],yt2[t],zt2[t],wt2[t],xt3[t],yt3[t],zt3[t],wt3[t],xt4[t],yt4[t],zt4[t],wt4[t]},{t,0,100},MaxSteps->100000];

A nyers ézrékenységek:

GraphicsArray[Table[Plot[lista2[[i,j]]/.mo/.refparam,{t,0,100},PlotLabel->{i,j}],{j,18},{i,4}]]

Periódusérzékenység :

S1[t]={{xa1[t],xa21[t],xa31[t],xa41[t],xa51[t],xa7[t],xa22[t],xa32[t],xa42[t],xa52[t],xl1[t],xl2[t],xj3[t],xj6[t],xt1[t],xt2[t],xt3[t],xt4[t]},{ya1[t],ya21[t],ya31[t],ya41[t],ya51[t],ya7[t],ya22[t],ya32[t],ya42[t],ya52[t],yl1[t],yl2[t],yj3[t],yj6[t],yt1[t],yt2[t],yt3[t],yt4[t]},{za1[t],za21[t],za31[t],za41[t],za51[t],za7[t],za22[t],za32[t],za42[t],za52[t],zl1[t],zl2[t],zj3[t],zj6[t],zt1[t],zt2[t],zt3[t],zt4[t]},{wa1[t],wa21[t],wa31[t],wa41[t],wa51[t],wa7[t],wa22[t],wa32[t],wa42[t],wa52[t],wl1[t],wl2[t],wj3[t],wj6[t],wt1[t],wt2[t],wt3[t],wt4[t]}}/.mo;

S1=S1[t]/.t->100;

SingularValueDecomposition[S1];

v={SingularValueDecomposition[S1]}[[1,3,1]]

sigma={SingularValueDecomposition[S1]}[[1,2,1,1]]

norma=Norm[modeleq[[All,2]]]/.refparam/.mo;

S1tau=tau*sigma/norma[[1]]/t*v;

Plot[Evaluate[S1tau[[1]]],{t,0,100}]

Normált periódusérzékenységek :

GraphicsArray[Table[Plot[S1tau[[i]]*i/tau/.mo/.refparam,{t,0,100},PlotLabel->{i}],{i,18}]]

Másik megoldás a periódusérzékenységek (Sτ) kiszámítására:

motnelkul=NDSolve[Join[modeleq,inipar,erzekeny,inierzekeny]/.refparam,{x,y,z,w,xa1,ya1,za1,wa1,xa21,ya21,za21,wa21,xa31,ya31,za31,wa31,xa41,ya41,za41,wa41,xa51,ya51,za51,wa51,xa7,ya7,za7,wa7,xa22,ya22,za22,wa22,xa32,ya32,za32,wa32,xa42,ya42,za42,wa42,xa52,ya52,za52,wa52,xl1,yl1,zl1,wl1,xl2,yl2,zl2,wl2,xj3,yj3,zj3,wj3,xj6,yj6,zj6,wj6,xt1,yt1,zt1,wt1,xt2,yt2,zt2,wt2,xt3,yt3,zt3,wt3,xt4,yt4,zt4,wt4},{t,0,100},MaxSteps->100000];

lista2tnelkul=lista2/.{xx_[t]->xx}

S1:=(lista2tnelkul[[1]]/.motnelkul)[[1]]

S1

szaml[t_]:=Through[S1[t]]-Through[S1[t+tau]]

nev[t_]:=modeleq[[1,2]]/.refparam/.nds

konst[i_][t_]:=szaml[t][[i]]/nev[t]

kozel[i_][delta_]:=Table[konst[i][t][[1]],{t,40,40+2*16.7346,delta}]

kozel[2][10]

int[i_][delta_,k_]:=Integrate[Interpolation[kozel[i][delta],InterpolationOrder->k][x],{x,40,40+2*16.7346}]

TableForm[Table[int[1][delta,3],{delta,0.0001,0.0010,0.0001}],TableHeadings->{Range[0.0001,0.0010,0.0001],None}]

stau[1]=860.5472057769115`/(2*16.7346)

Timing[TableForm[Table[int[2][delta,3],{delta,0.0001,0.0010,0.0001}],TableHeadings->{Range[0.0001,0.0010,0.0001],None}]]

stau[2]=-1088.604088875056`/(2*16.7346)

TableForm[Table[int[3][delta,3],{delta,0.0001,0.0010,0.0001}],TableHeadings->{Range[0.0001,0.0010,0.0001],None}]

stau[3]=4739.3370018618525`/(2*16.7346)

TableForm[Table[int[4][delta,3],{delta,0.0001,0.0010,0.0001}],TableHeadings->{Range[0.0001,0.0010,0.0001],None}]

stau[4]=-4818.222583936737`/(2*16.7346)

TableForm[Table[int[5][delta,3],{delta,0.0001,0.0010,0.0001}],TableHeadings->{Range[0.0001,0.0010,0.0001],None}]

stau[5]=-407.340334432147`/(2*16.7346)

TableForm[Table[int[6][delta,3],{delta,0.0001,0.0010,0.0001}],TableHeadings->{Range[0.0001,0.0010,0.0001],None}]

stau[6]=54.01621304815041`/(2*16.7346)

TableForm[Table[int[7][delta,3],{delta,0.0001,0.0010,0.0001}],TableHeadings->{Range[0.0001,0.0010,0.0001],None}]

stau[7]=-63.581673696734825`/(2*16.7346)

TableForm[Table[int[8][delta,3],{delta,0.0001,0.0010,0.0001}],TableHeadings->{Range[0.0001,0.0010,0.0001],None}]

stau[8]=158.2135493393998`/(2*16.7346)

TableForm[Table[int[9][delta,3],{delta,0.0001,0.0010,0.0001}],TableHeadings->{Range[0.0001,0.0010,0.0001],None}]

stau[9]=-164.4284020399711`/(2*16.7346)

TableForm[Table[int[10][delta,3],{delta,0.0001,0.0010,0.0001}],TableHeadings->{Range[0.0001,0.0010,0.0001],None}]

stau[10]=-892.5260084712627`/(2*16.7346)

TableForm[Table[int[11][delta,3],{delta,0.0001,0.0010,0.0001}],TableHeadings->{Range[0.0001,0.0010,0.0001],None}]

stau[11]=-0.15076454173289436`/(2*16.7346)

TableForm[Table[int[12][delta,3],{delta,0.0001,0.0010,0.0001}],TableHeadings->{Range[0.0001,0.0010,0.0001],None}]

stau[12]=13.620936820523848`/(2*16.7346)

TableForm[Table[int[13][delta,3],{delta,0.0001,0.0010,0.0001}],TableHeadings->{Range[0.0001,0.0010,0.0001],None}]

stau[13]=-10464.612817592768`/(2*16.7346)

TableForm[Table[int[14][delta,3],{delta,0.0001,0.0010,0.0001}],TableHeadings->{Range[0.0001,0.0010,0.0001],None}]

stau[14]=-5781.955652374687`/(2*16.7346)

TableForm[Table[int[15][delta,3],{delta,0.0001,0.0010,0.0001}],TableHeadings->{Range[0.0001,0.0010,0.0001],None}]

stau[15]=0.8200623536990118`/(2*16.7346)

TableForm[Table[int[16][delta,3],{delta,0.0001,0.0010,0.0001}],TableHeadings->{Range[0.0001,0.0010,0.0001],None}]

stau[16]=-401.63670381207476`/(2*16.7346)

TableForm[Table[int[17][delta,3],{delta,0.0001,0.0010,0.0001}],TableHeadings->{Range[0.0001,0.0010,0.0001],None}]

stau[17]=182.17670148414987`/(2*16.7346)

TableForm[Table[int[18][delta,3],{delta,0.0001,0.0010,0.0001}],TableHeadings->{Range[0.0001,0.0010,0.0001],None}]

stau[18]=469.03138694447495`/(2*16.7346)

A tisztított periódusérzékenység

GraphicsArray[Table[Plot[Evaluate[(lista2[[1,j]]+t/tau*stau[j]*(a1-(a21+a22*y[t])*x[t])/.refparam/.nds/.mo)],{t,0,100},PlotLabel->{1,j}],{j,18}]]

GraphicsArray[Table[Plot[Evaluate[(lista2[[2,j]]+t/tau*stau[j]*((((a31+a32*(t4-w[t]))*(t1-y[t]))/(j3+(t1-y[t])))-(((a41+a42*x[t])*y[t])/(j3+y[t])))/.refparam/.nds/.mo)],{t,0,100},PlotLabel->{2,j}],{j,18}]]

GraphicsArray[Table[Plot[Evaluate[(lista2[[3,j]]+t/tau*stau[j]*((((a51*(t4-w[t])+a52*t2)*(t3-z[t]-w[t]))/(j6+(t3-z[t]-w[t])))-((a7*x[t]*z[t])/(j6+z[t]))-l1*z[t]*(t4-w[t])+l2*w[t])/.refparam/.nds/.mo)],{t,0,100},PlotLabel->{3,j}],{j,18}]]

GraphicsArray[Table[Plot[Evaluate[(lista2[[4,j]]+t/tau*stau[j]*(l1*z[t]*(t4-w[t])-l2*w[t]-((a7*x[t]*w[t])/(j6+w[t])))/.refparam/.nds/.mo)],{t,0,100},PlotLabel->{4,j}],{j,18}]]

Tisztított normált periódusérzékenység

GraphicsArray[Table[Plot[Evaluate[((lista2[[1,j]]+t/tau*stau[j]*(a1-(a21+a22*y[t])*x[t]))/x[t]/.refparam/.nds/.mo)],{t,0,100},PlotLabel->{1,j}],{j,18}]]

GraphicsArray[Table[Plot[Evaluate[((lista2[[2,j]]+t/tau*stau[j]*((((a31+a32*(t4-w[t]))*(t1-y[t]))/(j3+(t1-y[t])))-(((a41+a42*x[t])*y[t])/(j3+y[t]))))/y[t]/.refparam/.nds/.mo)],{t,0,100},PlotLabel->{2,j}],{j,18}]]

GraphicsArray[Table[Plot[Evaluate[((lista2[[3,j]]+t/tau*stau[j]*((((a51*(t4-w[t])+a52*t2)*(t3-z[t]-w[t]))/(j6+(t3-z[t]-w[t])))-((a7*x[t]*z[t])/(j6+z[t]))-l1*z[t]*(t4-w[t])+l2*w[t]))/z[t]/.refparam/.nds/.mo)],{t,0,100},PlotLabel->{3,j}],{j,18}]]

GraphicsArray[Table[Plot[Evaluate[((lista2[[4,j]]+t/tau*stau[j]*(l1*z[t]*(t4-w[t])-l2*w[t]-((a7*x[t]*w[t])/(j6+w[t]))))/w[t]/.refparam/.nds/.mo)],{t,0,100},PlotLabel->{4,j}],{j,18}]]

Integrálközepek

vektorx[i_][t_]:=(lista2[[1,i]]+t/tau*stau[i]*(a1-(a21+a22*y[t])*x[t]))/x[t]/.refparam/.nds/.mo

vektory[i_][t_]:=(lista2[[2,i]]+t/tau*stau[i]*((((a31+a32*(t4-w[t]))*(t1-y[t]))/(j3+(t1-y[t])))-(((a41+a42*x[t])*y[t])/(j3+y[t]))))/y[t]/.refparam/.nds/.mo

vektorz[i_][t_]:=(lista2[[3,i]]+t/tau*stau[i]*((((a51*(t4-w[t])+a52*t2)*(t3-z[t]-w[t]))/(j6+(t3-z[t]-w[t])))-((a7*x[t]*z[t])/(j6+z[t]))-l1*z[t]*(t4-w[t])+l2*w[t]))/z[t]/.refparam/.nds/.mo

vektorw[i_][t_]:=(lista2[[4,i]]+t/tau*stau[i]*(l1*z[t]*(t4-w[t])-l2*w[t]-((a7*x[t]*w[t])/(j6+w[t]))))/w[t]/.refparam/.nds/.mo

kozelx[i_][delta_]:=Table[vektorx[i][t][[1]],{t,30,30+2*tau,delta}]

integralx[i_][delta_,k_]:=Integrate[Interpolation[kozelx[i][delta],InterpolationOrder->k][x],{x,30,30+2*tau}]

kozely[i_][delta_]:=Table[vektory[i][t][[1]],{t,30,30+2*tau,delta}]

integraly[i_][delta_,k_]:=Integrate[Interpolation[kozely[i][delta],InterpolationOrder->k][x],{x,30,30+2*tau}]

kozelz[i_][delta_]:=Table[vektorz[i][t][[1]],{t,30,30+2*tau,delta}]

integralz[i_][delta_,k_]:=Integrate[Interpolation[kozelz[i][delta],InterpolationOrder->k][x],{x,30,30+2*tau}]

kozelw[i_][delta_]:=Table[vektorw[i][t][[1]],{t,30,30+2*tau,delta}]

integralw[i_][delta_,k_]:=Integrate[Interpolation[kozelw[i][delta],InterpolationOrder->k][x],{x,30,30+2*tau}]

listax=Map[integralx[#][0.0001,3]&,Range[18]]

listay=Map[integraly[#][0.0001,3]&,Range[18]]

listaz=Map[integralz[#][0.0001,3]&,Range[18]]

listaw=Map[integralw[#][0.0001,3]&,Range[18]]

Az integrálközepek rendre :

lista=Join[listax,listay,listaz,listaw,{},{},{},{},{},{}]/(2tau)

Az érzékenységi vektor hossza :

Norm[lista]