# Solution of a system of differential equations with a correction factor varying over time

Dear Mathematica users,

I am trying to solve a system of equations of first order and time depedent. The system has 5 equations, one for every variable which varies over time. At the moment I am using NDSolve and every thing works fine. However, now my system got complicated. I do have to correct it by empliying a correction factor when my temperature increases above certain value. The temperature is one of the variables which depends of the time and the correction factor affects only to some parameters of certain equations. This correction factor changes from 1 to 10E-45. My fundamental problem is that I do not how to introduce this correction meanwhile the calculation is being made. I do not know whether I should use conditional and in the way a must use them. The solution would be some subroutine in which the system is solve for every time step and termperature is compare with the critical value. Thus the NDSolve can either continue and keep correction=1 or change the value to 10E-45.

Does anybody have any idea about how to solve this programming problem?

Thank you very much in advance. Regards.

Dear StackExchange Mathematica comunity, First of all, sorry for not introducing my system of equations in the question with the value of parameters. I thought it was going to be a lot. Secondly, sorry for introducing it in an answer and in a bad way. I do have my equations copied from Math as Plain Text and pasted in a Word file. From there one can copy them and write paste it in Math. It does work. So I did paste it again here in the same format.

Thank you very much again in advance (I do think that WhenEvent will work out, but please tell more possibilities, so I could learn about this fantastic program)

    (*Laser amplitude calculation*)
fluence=20*10^-3 *10000/1(*/J m-2*);
lspot=0.0001(*Laser spot diameter/ m*);
larea=(lspot*lspot/4)*Pi;
energy=fluence*(lspot*lspot/4)*Pi;(*/J*)
pulseduration=5*10^-9;(*/s*)
peakintensity=energy/pulseduration ;
intensity=fluence/pulseduration;
lim1=0;
lim2=50*10^-9;
lim3=100*10^-9;
μ1=8*10^-9;σ1=(4*10^-9)/2.3528;y5=0;
f[y_,r1_]:=y5+r1*(E^(( -(y-μ1)^2/(2 σ1^2)))) /(σ1 Sqrt[2 π]);
res=\!$$\*SubsuperscriptBox[\(∫$$, $$lim1$$, $$lim2$$]$$\((y5 + r1* \*FractionBox[\( \*SuperscriptBox[\(E$$,
FractionBox[$$\(\$$$$- \*SuperscriptBox[\((y - μ1)$$, $$2$$]\)\), $$2\ \*SuperscriptBox[\(σ1$$, $$2$$]\)]]$$\$$\), $$σ1 \*SqrtBox[\(2\ π$$]\)])\) \[DifferentialD]y\)\);
RES=NSolve[res==intensity,r1]
a2=r1/.RES[[1]];
a3=a2/(pulseduration*(lspot*lspot/4)*Pi);
a4=fluence/pulseduration;

(*Gaussian laser pulse*)
u[t_]:=y5+a2*(E^(( -(t-μ1)^2/(2 σ1^2)))) /(σ1 Sqrt[2 π]);
a5=1.602177*10^2;
(*Voltage function*)
del=10*10^-9;v=5000;τ3=1*10^-9;
z[y_]:=v-E^((- y+del)/τ3);
j[y_]:=1/(1+E^((- y+del)/τ3))

(*Parameters*)
p1=8.05*(1.602177*10^-19)/1;
p2=7.82*(1.602177*10^-19)/1(*7.2*);
n1=3.466*(1.602177*10^-19)/1;
n2=2*3.6789*(1.602177*10^-19)/1;
temp=298;
σ01=(9.9*10^-18)*1/10000;
σ10=(2*10^-18)*1/10000;
σ1n=(2*10^-18)*1/10000;
τ1=(1*10^-9-3*10^-9)/(1+E^((t-4.44*10^-9)/10^-15))+3*10^-9;
τ2=3*10^-9;
ϕs1=(1*10^-9)/(3*10^-9);
f=4.5;
density=1.44*1000000/1000;
kn1=7.5*10^10;
k11=7*10^9;
k1n=1*10^11;
k20=1.5*10^11;
fcluster=0.85;
dip=6;(*Range of pooling process*)
h=6.626*10^-34;
υ=(3*10^8)/(337*10^-9)(*Laser wavelength / m*);
time=1*10^-9;
b=1.38066*10^-23;
pansion=4;
conc=0.000005/1*6.23*10^23*10^-6
etial=temp*b*f *conc(*Internal energy*)

(*k thermic*)
k[m_]:=9*10^15*Exp[(n1-p2)/(b*m)];

na=1 (*6.23*10^23*);
matrixspot=0.0015^2*Pi;
areaspot=(lspot*lspot/4)*Pi;
ini=3.46*10^15(*conc*areaspot/matrixspot//(2.55744*10^5)/(360*10^2*10^-27) 0.05*10^-6*6.23*10^23*areaspot/matrixspot*)

(*Expansion Correction*)
c=6.982;
b1=1.489;
b2=3.036;
c1=1.5;
c2=0.4282;
lspot=0.0001;
ddistance=10^-7;
equis=ddistance/lspot;
mfactor=c-b1*Exp[-c1*equis]-b2*Exp[-c2*equis];
γ=1.05;
tau=(4*10^-9)/2.3528;

correction=(1+γ/(γ-1)*mfactor^2)^(-γ/(γ-1));
corrt[t_]:=Piecewise[{{1,t<4.43*10^-9},{correction,t>=4.43*10^-9}}]

• Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. – bbgodfrey May 6 '15 at 21:05
• Please include a simple version (say, one equation) of your problem, including the correction factor, so that readers can make concrete suggestions. – bbgodfrey May 6 '15 at 21:07
• If you are using MMA of v.9 or higher, I can suggest looking up WhenEvent in the documentation. – LLlAMnYP May 6 '15 at 23:07

It is unclear where the correction factor is to be applied in the OP's code, or at least it is on a cursory reading of the code. The problem in the title is as easy as defining a function and including it as a coefficient in the differential equation:

correctionfactor[t_] := ...
NDSolve[{x'[t] == correctionfactor[t] x[t]...


The general problem as described in the text of the question, in which the coefficient is discontinuous, needs special treatment for NDSolve to efficiently handle the discontinuities; a solution can be implemented with either Piecewise or WhenEvent.

Piecewise changes the value of correctionfactor whenever the temperature T crosses the criticaltemp. Here the ODE defines temp[t] as a monotonic function, so this happens at most once.

Clear[correctionfactor];
criticaltemp = 1;
correctionfactor[T_] := Piecewise[{{1, T < criticaltemp}}, 10^-1];
ode = {temp'[t] == correctionfactor[temp[t]] temp[t], temp[0] == 0.1};
{sol} = NDSolve[ode, temp, {t, 0, 10}];

Plot[temp[t] /. sol // Evaluate, {t, 0, 10}]


This WhenEvent changes the correctionfactor only when the temperature increases past the criticaltemp. Since temp[t] is monotonic, this is equivalent to the first method.

Clear[correctionfactor];
criticaltemp = 1;
ode = {temp'[t] == correctionfactor[t] temp[t], temp[0] == 0.1};
{sol} = NDSolve[{ode,
correctionfactor[0] == 1,
WhenEvent[temp[t] > criticaltemp, correctionfactor[t] -> 10^-1]},
temp, {t, 0, 10},
DiscreteVariables -> {correctionfactor}];

Plot[temp[t] /. sol // Evaluate, {t, 0, 10}]


Below the ODE leads to an oscillating temp[t]. With Piecewise we see the discontinuity in the derivative every time the temperature crosses (whether increasing or decreasing) the criticaltemp.

Clear[correctionfactor];
criticaltemp = 1;
correctionfactor[t_] := Piecewise[{{1, t < criticaltemp}}, 10^-45];
ode = {temp'[t] == -correctionfactor[temp[t]] temp[t] + 2 Cos[t], temp[0] == 0.2};
{sol} = NDSolve[ode, temp, {t, 0, 25}];

Plot[temp[t] /. sol // Evaluate, Evaluate@Flatten[{t, temp["Domain"] /. sol}]]


But the single WhenEvent changes the correctionfactor to 10^-45 and never changes it back:

Clear[correctionfactor];
criticaltemp = 1;
ode = {temp'[t] == -correctionfactor[t] temp[t] + 2 Cos[t], temp[0] == 0.1};
{sol} = NDSolve[{ode,
correctionfactor[0] == 1,
WhenEvent[temp[t] > criticaltemp, correctionfactor[t] -> 10^-45]},
temp, {t, 0, 25},
DiscreteVariables -> {correctionfactor}];

Plot[temp[t] /. sol // Evaluate, Evaluate@Flatten[{t, temp["Domain"] /. sol}]]


To get the same behavior as Piecewise, we need a second WhenEvent.

Clear[correctionfactor];
criticaltemp = 1;
ode = {temp'[t] == -correctionfactor[t] temp[t] + 2 Cos[t], temp[0] == 0.1};
{sol} = NDSolve[{ode,
correctionfactor[0] == 1,
WhenEvent[temp[t] > criticaltemp, correctionfactor[t] -> 10^-45],
WhenEvent[temp[t] < criticaltemp, correctionfactor[t] -> 1]},
temp, {t, 0, 25},
DiscreteVariables -> {correctionfactor}];

Plot[temp[t] /. sol // Evaluate, Evaluate@Flatten[{t, temp["Domain"] /. sol}]]


Which of the last two WhenEvent models is appropriate depends entirely on the problem being modeled. WhenEvent has options that makes building a robust model somewhat easier in some cases than Piecewise.

Caveat: Note that the swing from 1 to 10^-45 is quite large at ordinary MachinePrecision, which has slightly less than 16 digits of precision. It might lead to problems with the numerics. See also \$MachineEpsilon.

• Dear Michael E2, – RKdifferential May 10 '15 at 11:50
• First of all, thank´s a lot for your comment. I will try this way write now. Besides, I did change the equations so, you can see where the correction must be applied. Only to some constants (K20,K11,K10).Secondly I would like to ask you for something more. I should repeat the calculation several time for different ddistance values (see Expansion correction in the text) and then I should sum up all the valued i go from every integration. How can I do this? In other programming codes one can used Do, but in Mathematica this does not work out. What do you recommend me? Thank you in advance! – RKdifferential May 10 '15 at 12:03
• Dear Michael E2, I have implemented piecewise for my correction factor and for two other parameters that change over time. Now my problem is another. If I set a StartingStepSize lower than te time (10E-9) where my T=450, the program stops and returns: – RKdifferential May 10 '15 at 14:18
• NDSolve::ndsz: At t == 1.034911650603046`*^-8, step size is effectively zero; singularity or stiff system suspected. >> – RKdifferential May 10 '15 at 14:18
• Do you Know how can I solve this? – RKdifferential May 10 '15 at 14:18