0
$\begingroup$

I would like to solve a set of differential equations but vary one parameter, KieP, between 900 and 1100.

The solution for a given parameter has many peaks. I would then like to interpolate between the peaks to create a smooth envelope.

Finally, I would like to integrate the envelope and plot this value against the parameter, KieP.

I have written some code which works for a single parameter:

(*Defining Parameters*)
Inj[t_] := Piecewise[{{0, t < 0}, {1, 0 < t < 10}, {0, t > 10}}]; 
InjCoeff1 = 1;
InjCoeff2 = 1; 
T1PE = 40; 
T1PI = 60; 
T2P = 1/60;
T1LE = 40; 
T1LI = 60; 
T2L = 1/60;
Klp = 0.004;
Kpl = 0.004;
KieP = 1000; 
KeiP = 1;  
KieL = 1;  
KeiL = 0.1;
\[Epsilon] = 0.0001 ;
\[Alpha] = 5;          
\[Theta]z[t_] := \!\(
 \*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(500\)]\(
  \*FractionBox[
   SuperscriptBox[\(E\), 
    FractionBox[\(-\ 
      \*SuperscriptBox[\((t - 
         i)\), \(2\)]\), \(4  \[Epsilon]\)]], \(2 
    \*SqrtBox[\(\[Epsilon]*\[Pi]\)]\)] \((Cos[\[Alpha]\ Degree] - 
    1)\)\)\);

\[Theta]xy[t_] := \!\(
 \*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(500\)]\(
  \*FractionBox[
   SuperscriptBox[\(E\), 
    FractionBox[\(-\ 
      \*SuperscriptBox[\((t - 
         i)\), \(2\)]\), \(4  \[Epsilon]\)]], \(2 
    \*SqrtBox[\(\[Epsilon]*\[Pi]\)]\)] \((Sin[\[Alpha]\ \
Degree])\)\)\);

(*Equations*)
sol = NDSolve[{Pez'[t] == 
InjCoeff1*Inj[t] - 1/T1PE*Pez[t] + KieP*Piz[t] - 
 KeiP*Pez[t] + \[Theta]z[t]*Pez[t], 
Piz'[t] ==  -(1/T1PI)*Piz[t] + Klp*Liz[t] - Kpl*Piz[t] + 
 KeiP*Pez[t] - KieP*Piz[t] + \[Theta]z[t]*Piz[t],
Lez'[t] ==  -(1/T1LE)*Lez[t] + KieL*Liz[t] - 
 KeiL*Lez[t] + \[Theta]z[t]*Lez[t],
Liz'[t] ==  -(1/T1LI)*Liz[t] + Kpl*Piz[t] + KeiL*Lez[t] - 
 Klp*Liz[t] - KieL*Liz[t] + \[Theta]z[t]*Liz[t], 
Pexy'[t] == -(1/T2P)*Pexy[t] + KieP*Pixy[t] - 
 KeiP*Pexy[t] + \[Theta]xy[t]*Pez[t - 5*Sqrt[2*\[Epsilon]]],
Pixy'[t] ==  -(1/T2P)*Pixy[t] + Klp*Lixy[t] - Kpl*Pixy[t] + 
 KeiP*Pexy[t] - 
 KieP*Pixy[t] + \[Theta]xy[t]*Piz[t - 5*Sqrt[2*\[Epsilon]]], 
Lexy'[t] == -(1/T2L)*Lexy[t] + KieL*Lixy[t] - 
 KeiL*Lexy[t] + \[Theta]xy[t]*Lez[t - 5*Sqrt[2*\[Epsilon]]], 
Lixy'[t] == -(1/T2L)*Lixy[t] + Kpl*Pixy[t] - Klp*Lixy[t] + 
 KeiL*Lexy[t] - 
 KieL*Lixy[t] + \[Theta]xy[t]*Liz[t - 5*Sqrt[2*\[Epsilon]]], 
Pet'[t] == 
InjCoeff2*Inj[t] + 1/T1PE*Pez[t] + 1/T2P*Pexy[t] + KieP*Pit[t] - 
 KeiP*Pet[t], 
Pit'[t] == 
1/T1PI*Piz[t] + 1/T2P*Pixy[t] + Klp*Lit[t] - Kpl*Pit[t] + 
 KeiP*Pet[t] - KieP*Pit[t],
Let'[t] ==  
1/T1LE*Lez[t] + 1/T2L*Lexy[t] + KieL*Lit[t] - KeiL*Let[t],
Lit'[t] == 
1/T1LI*Liz[t] + 1/T2L*Lixy[t] + Kpl*Pit[t] - Klp*Lit[t] + 
 KeiL*Let[t] - KieL*Lit[t],

Pez[0] == 0, Piz[0] == 0, Lez[0] == 0, Liz[0] == 0, Pexy[0] == 0 , 
Pixy[0] == 0, Lexy[0] == 0, Lixy[0] == 0, Pet[0] == 0.1, 
Pit[0] == 0.15, Let[0] == 1, Lit[0] == 3.5}, 
{Pez , Piz,   Lez,   Liz,   Pexy, Pixy, Lexy, Lixy, Pet,  Pit,   
Let,   Lit , Ltotxy, Ptotxy  }, {t, 0, 250}]


(*Defining New varaiables*)
Ptotxy[t_] := Evaluate[Pexy[t] /. sol] + Evaluate[Pixy[t] /. sol]
Ltotxy[t_] := Evaluate[Lexy[t] /. sol] + Evaluate[Lixy[t] /. sol]
(*Taking the derivative*)
A[t_] := Ptotxy'[t]
B[t_] := Ltotxy'[t]

(*Finding the peaks of the graphs*)
maximaP = Flatten[FindRoot[A[t], {t, Table[i, {i, 249}]}][[1, 2]]];
maximaL = Flatten[FindRoot[B[t], {t, Table[i, {i, 249}]}][[1, 2]]];

(*Creating an series of coordiates*)
datapointsP = Riffle[maximaP, Flatten[Ptotxy[maximaP]]]~Partition~2;
datapointsL = Riffle[maximaL, Flatten[Ltotxy[maximaL]]]~Partition~2;

(*Interpolation*)
datapointsPInterp = Interpolation[datapointsP];
datapointsLInterp = Interpolation[datapointsL];

(*Numerical Integration*)
AUCP = NIntegrate[datapointsPInterp[t], {t, 0, 249}];
AUCL = NIntegrate[datapointsLInterp[t], {t, 0, 249}];
AUCRatio = AUCL/AUCP

I have tried using ParametricNDSolve but have got totally lost with the syntax.

I thought a simpler way may be to use the map function. Here is what I have done so far:

I have left the

(*Defining Parameters*)

section the same except change KieP to

KiePrange = Range[900, 1100, 50];

Then changing KieP to # in NDSolve:

sol = NDSolve[{Pez'[t] == 
   InjCoeff1*Inj[t] - 1/T1PE*Pez[t] + #*Piz[t] - 
    KeiP*Pez[t] + \[Theta]z[t]*Pez[t], 
  Piz'[t] ==  -(1/T1PI)*Piz[t] + Klp*Liz[t] - Kpl*Piz[t] + 
    KeiP*Pez[t] - #*Piz[t] + \[Theta]z[t]*Piz[t],
  Lez'[t] ==  -(1/T1LE)*Lez[t] + KieL*Liz[t] - 
    KeiL*Lez[t] + \[Theta]z[t]*Lez[t],
  Liz'[t] ==  -(1/T1LI)*Liz[t] + Kpl*Piz[t] + KeiL*Lez[t] - 
    Klp*Liz[t] - KieL*Liz[t] + \[Theta]z[t]*Liz[t], 
  Pexy'[t] == -(1/T2P)*Pexy[t] + #*Pixy[t] - 
    KeiP*Pexy[t] + \[Theta]xy[t]*Pez[t - 5*Sqrt[2*\[Epsilon]]],
  Pixy'[t] ==  -(1/T2P)*Pixy[t] + Klp*Lixy[t] - Kpl*Pixy[t] + 
    KeiP*Pexy[t] - #*Pixy[t] + \[Theta]xy[t]*
     Piz[t - 5*Sqrt[2*\[Epsilon]]], 
  Lexy'[t] == -(1/T2L)*Lexy[t] + KieL*Lixy[t] - 
    KeiL*Lexy[t] + \[Theta]xy[t]*Lez[t - 5*Sqrt[2*\[Epsilon]]], 
  Lixy'[t] == -(1/T2L)*Lixy[t] + Kpl*Pixy[t] - Klp*Lixy[t] + 
    KeiL*Lexy[t] - 
    KieL*Lixy[t] + \[Theta]xy[t]*Liz[t - 5*Sqrt[2*\[Epsilon]]], 
  Pet'[t] == 
   InjCoeff2*Inj[t] + 1/T1PE*Pez[t] + 1/T2P*Pexy[t] + #*Pit[t] - 
    KeiP*Pet[t], 
  Pit'[t] == 
   1/T1PI*Piz[t] + 1/T2P*Pixy[t] + Klp*Lit[t] - Kpl*Pit[t] + 
    KeiP*Pet[t] - #*Pit[t],
  Let'[t] ==  
   1/T1LE*Lez[t] + 1/T2L*Lexy[t] + KieL*Lit[t] - KeiL*Let[t],
  Lit'[t] == 
   1/T1LI*Liz[t] + 1/T2L*Lixy[t] + Kpl*Pit[t] - Klp*Lit[t] + 
    KeiL*Let[t] - KieL*Lit[t],

  Pez[0] == 0, Piz[0] == 0, Lez[0] == 0, Liz[0] == 0, 
  Pexy[0] == 0 , Pixy[0] == 0, Lexy[0] == 0, Lixy[0] == 0, 
  Pet[0] == 0.1, Pit[0] == 0.15, Let[0] == 1, Lit[0] == 3.5}, 
 {Pez , Piz,   Lez,   Liz,   Pexy, Pixy, Lexy, Lixy, Pet,  Pit,   
  Let,   Lit , Ltotxy, Ptotxy  }, {t, 0, 250}] & /@ KiePrange;

Ptotxy[t_] := Evaluate[Pexy[t] /. #] + Evaluate[Pixy[t] /. #] & /@ sol;
Ltotxy[t_] := Evaluate[Lexy[t] /. #] + Evaluate[Lixy[t] /. #] & /@ sol;
(*Taking the derivative*)
A[t_] := Ptotxy'[t];
B[t_] := Ltotxy'[t];

From here, I am unsure how to adapt my code to get it to solve for the different values of KieP with the end goal to plot AUCRatio vs KieP. Thanks in advance for any help!

$\endgroup$
  • $\begingroup$ If you only want to vary one parameter in a numerically solution ParametricNDSolve could solve your problem $\endgroup$ – Ulrich Neumann Sep 11 '18 at 15:21
  • $\begingroup$ @UIrich have tried, but I have found the syntax very confusing when doing all of the data manipulations after ParametricNDSolve. $\endgroup$ – Dave Bassito Sep 11 '18 at 15:38
  • 1
    $\begingroup$ Is this simpler and less confusing? Table[(sol = NDSolve[...];Ptotxy[t_] :=...;AUCRatio = AUCL/AUCP; {AUCRatio, KieP}), {KieP, 900, 1100, 50}] That puts your whole solution process inside () inside Table so all of that code will be interpreted for every step. That gives you a list of pairs of the two variables you are want. Then you can ListPlot that list. When I do that I get the list of five pairs of values along with a variety of warnings that are all from your NDSolve that seem to have nothing to do with how you then post-process your solution. $\endgroup$ – Bill Sep 11 '18 at 16:04
  • $\begingroup$ Thanks @Bill, this has really helped! $\endgroup$ – Dave Bassito Sep 11 '18 at 17:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.