1
$\begingroup$

Main Question:

  1. Why does NDSolve fail with Power::infy in Versions 14.3.0.0, 14.2.1.0, and 14.2.0.0, yet succeeds in Version 14.0.0.0?
  2. More troubling, why does NDSolve fail with Power::infy in Version 14.0.0.0 when I change the names of the variables? This suggests that the distinguishing factor is the order in which NDSolve process the variables.
  3. Is the naming of the variables the only mechanism for influencing the order in which NDSolve process the variables? If so, what rules dictate the processing order? If not, what other mechanisms are available?
  4. Is there a better way to avoid these errors?

Background:

  • In Versions 14.3.0.0, 14.2.1.0, and 14.2.0.0, NDSolve terminates with Power::infy and the message "Encountered non-numerical value for a derivative at t == 1.`60."
  • In Version 14.0.0.0, NDSolve successfully solves the problem
  • When I change the names of the variables ([DoubledGamma] ->obj and [Lambda]->w), even in Version 14.0.0.0, NDSolve terminates Power::infy and the message "Encountered non-numerical value for a derivative at t == 1.`60."

Details: Here's the smallest and simplest example I can construct. I am trying to maintain the conditions:

k[1][t]-k[9][t]+Integral from 0 to s[1,9][t] of (f[x, 1, t] - f[x, 9, t])dx = K[1,9]+Epsilon*t,

k[5][t]-k[9][t]+Integral from 0 to s[5,9][t] of (f[x, 5, t] - f[x, 9, t])dx = K[5,9]+Epsilon*t,

g[0,9,t]=-1,

g[0,10,t]=-1,

g[1,1,t]=2,

g[1,5,t]=2,

k[9][t]-k[5][t]+Integral from 0 to s[9,5][t] of (f[x, 9, t] - f[x, 5, t])dx = K[9,5]+Epsilon*t,

k[10][t]-k[5][t]+Integral from 0 to s[10,5][t] of (f[x, 10, t] - f[x, 5, t])dx = K[10,5]+Epsilon*t,

k[5][t]-k[10][t]+Integral from 0 to z[10,5][t] of (f[x, 5, t] - f[x, 10, t])dx = K[5,10]+Epsilon*t,

while reducing t from the initial value 1 to 0.

Satisfying these conditions ensures \[DoubledGamma] is as small as possible for the given value of Epsilon. They are satisfied at t=1. We want to find \[DoubledGamma][0], the smallest value of \[DoubledGamma] with Epsilon = 0.

The functions g are defined by:

g[x_,i_,t_]:=1+1/(2 u[i]^2)- x/u[i]+\[DoubledGamma][t]/u[i]+E^(-2 x u[i]) c[i][t];

for i=1,5,9,10 and

rates = {u[1]->-5, u[2]->-4, u[3]->-3, u[4]->-2, u[5]->-1, u[9]->4, u[10]->5};

and the functions f are defined by

f[x_,1,t_] := Piecewise[{{g[x, 10, t], 0 <= x <= s[1, 10][t]}, {g[x, 1, t], s[1, 10][t] <= x <= 1}}, 0];

f[x,5,t] is the convex combination:

f[x_,5,t] := Piecewise[{{g[x, 9, t], 0 <= x <= s[5, 9][t]}, {g[x, 5, t], s[5, 9][t] <= x <= 1}}, 0]*\[Lambda][{5, {9, 5}}][t] + Piecewise[{{g[x, 10, t], 0 <= x <= s[5, 10][t]}, {g[x, 5, t], s[5, 10][t] <= x <= 1}}, 0]*\[Lambda][{5, {10, 5}}][t];

So, \[Lambda][{5,{9,5}}][t] and \[Lambda][{5,{10,5}}][t] are non-negative and add to 1. Finally, f[x,9,t] and f[x,10,t] are defined by:

f[x_,i_/;i>8,t_] := Piecewise[{{g[x, i, t], 0 <= x <= s[i, 1][t]}, {g[x, 1, t], s[i, 1][t] <= x <= 1}}, 0]

The points s[i,j][t] satisfy g[s[i,j][t],i,t]=g[s[i,j][t],j,t]. The two point s[i,j][t] and s[j,i][t] are different points of intersection. The point z[5,10][t] satisfies f[z[5,10][t],5,t]=f[z[5,10][t],10,t]. In this case it coincides with s[9,10][t]=0, but in general that need not be the case.

These points of intersection motivate the use of NDSolve, which must compute the system variable \[DoubledGamma][t], the shape scalars c[i][t], the distance scalars k[i][t] and the weights \[Lambda][{5,{9,5}}][t] and \[Lambda][{5,{10,5}}][t] as well as the points s[i,j][t] and z[5,10][t] as t decreases.

NDSolve generally works well. But in this case encounters issues finding initial expressions for the derivatives. High precision (60) is required here.

NDSolve must solve for the main variables: \[DoubledGamma], c, k, the weights \[Lambda] and the points of intersection s and z:

mainVars = {\[DoubledGamma],c[1],c[5],c[9],c[10],k[1],k[5],k[9],k[10],\[Lambda][{1,{10,1}}],\[Lambda][{5,{9,5}}],\[Lambda][{9,{9,1}}],\[Lambda][{10,{10,1}}],\[Lambda][{5,{10,5}}]}; sVars = {s[5,1],s[1,5],s[9,1],s[1,9],s[9,5],s[5,9],s[10,1],s[1,10],s[10,5],s[5,10],s[10,9],s[9,10]}; zVars = {z[5,10]};

The initial values are:

initvals = {\[DoubledGamma][1]==2.05447311893110626081720852402808538099901372370956350616897411079436318649221`59.21203534568774,\[Lambda][{1,{10,1}}][1]==1.`59.52287874528034,\[Lambda][{5,{9,5}}][1]==1.`59.52287874528034,\[Lambda][{9,{9,1}}][1]==1.`59.52287874528034,\[Lambda][{10,{10,1}}][1]==1.`59.52287874528034,\[Lambda][{5,{10,5}}][1]==0,c[1][1]==0.00005406653227441526643868237352926323466039253045343914380588423604120768158`59.38908429919841,c[5][1]==0.21037505983424199245791900612097248467182138429380183947099074315029886791301`59.14592587330862,c[9][1]==-2.54486827973277656520430213100702134524975343092739087654224352823976711804799`59.43971671306184,c[10][1]==-2.43089462378622125216344170480561707619980274474191270123379482259181369444374`59.45219230677731,k[1][1]==0.5408776874057795130454721092380356600764207541942596435546875`59.52287874528034,k[5][1]==0.55465138667532199656332331136023086604109198780433921296205222358427469602437`59.510756214539754,k[9][1]==0.27715144945091994400932213001761150975529731837187564925111034467168456281922`58.9691397433357,k[10][1]==0.27460703309891994424949125726005461398097382512961185530579850730766021922944`59.00514896990515,s[5,1][1]==1,s[1,5][1]==0.60349555274901874397652223337236012112694574140435952764373361673556150988203`59.308361932635,s[9,1][1]==0.91738402987730429611143918835794537392089729246177868444942596314298054187187`59.46145928246406,s[1,9][1]==0.13337426812926725479441414451141974687176300977327092635185534286252287751542`59.249336965538,s[9,5][1]==0.81592390857510720258798682321748712700851541037656172491318302135643499946369`59.36413379244892,s[5,9][1]==0.04204597590699506909254996282792006963593968045792112547777783770548567614268`58.85860660094471,s[10,1][1]==0.90478712690090498019741225952927798672556918947223242437794926149680359450133`59.46018924130699,s[1,10][1]==0.11419363182898146466115313366387218129562463476248516091019311833200443800258`59.257500221171576,s[10,5][1]==0.79465702812305202254114584951192129601351644032525243179532741600110321054715`59.360946380290144,s[5,10][1]==0.03469582551727117232910698890610149243752661355165802816814932543205104711335`58.862209858860176,s[10,9][1]==0.3193782270321289973211897540302884699083663997202095952809160860280949029792`59.35316694806237,s[9,10][1]==0,z[5,10][1]==0};

The differential equations to maintain the conditions above are:

mainEqns = {Derivative[1][k[1]][t]==0,Derivative[1][\[Lambda][{1,{10,1}}]][t]==0,Derivative[1][\[Lambda][{5,{9,5}}]][t]+Derivative[1][\[Lambda][{5,{10,5}}]][t]==0,Derivative[1][\[Lambda][{9,{9,1}}]][t]==0,Derivative[1][\[Lambda][{10,{10,1}}]][t]==0,-\[Lambda][{9,{9,1}}][t] (1/4 Min[s[1,9][t],s[9,1][t]] Derivative[1][\[DoubledGamma]][t]+1/8 Derivative[1][c[9]][t]-1/8 E^(-8 Min[s[1,9][t],s[9,1][t]]) Derivative[1][c[9]][t])+\[Lambda][{1,{10,1}}][t] (-(1/5) Min[1,s[1,9][t]] Derivative[1][\[DoubledGamma]][t]+1/5 Min[s[1,9][t],s[1,10][t]] Derivative[1][\[DoubledGamma]][t]+1/5 s[1,10][t] Derivative[1][\[DoubledGamma]][t]+1/10 E^(10 Min[1,s[1,9][t]]) Derivative[1][c[1]][t]-1/10 E^(10 s[1,10][t]) Derivative[1][c[1]][t]+1/10 Derivative[1][c[10]][t]-1/10 E^(-10 Min[s[1,9][t],s[1,10][t]]) Derivative[1][c[10]][t])+Derivative[1][k[1]][t]-Derivative[1][k[9]][t]+(51/50 Min[1,s[1,9][t]]+1/10 Min[1,s[1,9][t]]^2+51/50 Min[s[1,9][t],s[1,10][t]]-1/10 Min[s[1,9][t],s[1,10][t]]^2-1/5 Min[1,s[1,9][t]] \[DoubledGamma][t]+1/5 Min[s[1,9][t],s[1,10][t]] \[DoubledGamma][t]+1/10 E^(10 Min[1,s[1,9][t]]) c[1][t]-1/10 E^(10 s[1,10][t]) c[1][t]+1/10 c[10][t]-1/10 E^(-10 Min[s[1,9][t],s[1,10][t]]) c[10][t]-51/50 s[1,10][t]+1/5 \[DoubledGamma][t] s[1,10][t]-1/10 s[1,10][t]^2) Derivative[1][\[Lambda][{1,{10,1}}]][t]-(33/32 Min[s[1,9][t],s[9,1][t]]-1/8 Min[s[1,9][t],s[9,1][t]]^2+1/4 Min[s[1,9][t],s[9,1][t]] \[DoubledGamma][t]+1/8 c[9][t]-1/8 E^(-8 Min[s[1,9][t],s[9,1][t]]) c[9][t]) Derivative[1][\[Lambda][{9,{9,1}}]][t]==0.276281253576402052622079925762227839726037377161803163320316878162483075729117663799782079680408784210336126797460477138479369335828285`121.62355863066968,-\[Lambda][{9,{9,1}}][t] (1/4 Min[s[5,9][t],s[9,1][t]] Derivative[1][\[DoubledGamma]][t]+1/8 Derivative[1][c[9]][t]-1/8 E^(-8 Min[s[5,9][t],s[9,1][t]]) Derivative[1][c[9]][t])+\[Lambda][{5,{9,5}}][t] (-Min[1,s[5,9][t]] Derivative[1][\[DoubledGamma]][t]+5/4 s[5,9][t] Derivative[1][\[DoubledGamma]][t]+1/2 E^(2 Min[1,s[5,9][t]]) Derivative[1][c[5]][t]-1/2 E^(2 s[5,9][t]) Derivative[1][c[5]][t]+1/8 Derivative[1][c[9]][t]-1/8 E^(-8 s[5,9][t]) Derivative[1][c[9]][t])+\[Lambda][{5,{10,5}}][t] (-Min[1,s[5,9][t]] Derivative[1][\[DoubledGamma]][t]+1/5 Min[s[5,9][t],s[5,10][t]] Derivative[1][\[DoubledGamma]][t]+s[5,10][t] Derivative[1][\[DoubledGamma]][t]+1/2 E^(2 Min[1,s[5,9][t]]) Derivative[1][c[5]][t]-1/2 E^(2 s[5,10][t]) Derivative[1][c[5]][t]+1/10 Derivative[1][c[10]][t]-1/10 E^(-10 Min[s[5,9][t],s[5,10][t]]) Derivative[1][c[10]][t])+Derivative[1][k[5]][t]-Derivative[1][k[9]][t]+(3/2 Min[1,s[5,9][t]]+1/2 Min[1,s[5,9][t]]^2-Min[1,s[5,9][t]] \[DoubledGamma][t]+1/2 E^(2 Min[1,s[5,9][t]]) c[5][t]-1/2 E^(2 s[5,9][t]) c[5][t]+1/8 c[9][t]-1/8 E^(-8 s[5,9][t]) c[9][t]-15/32 s[5,9][t]+5/4 \[DoubledGamma][t] s[5,9][t]-5/8 s[5,9][t]^2) Derivative[1][\[Lambda][{5,{9,5}}]][t]+(3/2 Min[1,s[5,9][t]]+1/2 Min[1,s[5,9][t]]^2+51/50 Min[s[5,9][t],s[5,10][t]]-1/10 Min[s[5,9][t],s[5,10][t]]^2-Min[1,s[5,9][t]] \[DoubledGamma][t]+1/5 Min[s[5,9][t],s[5,10][t]] \[DoubledGamma][t]+1/2 E^(2 Min[1,s[5,9][t]]) c[5][t]-1/2 E^(2 s[5,10][t]) c[5][t]+1/10 c[10][t]-1/10 E^(-10 Min[s[5,9][t],s[5,10][t]]) c[10][t]-3/2 s[5,10][t]+\[DoubledGamma][t] s[5,10][t]-1/2 s[5,10][t]^2) Derivative[1][\[Lambda][{5,{10,5}}]][t]-(33/32 Min[s[5,9][t],s[9,1][t]]-1/8 Min[s[5,9][t],s[9,1][t]]^2+1/4 Min[s[5,9][t],s[9,1][t]] \[DoubledGamma][t]+1/8 c[9][t]-1/8 E^(-8 Min[s[5,9][t],s[9,1][t]]) c[9][t]) Derivative[1][\[Lambda][{9,{9,1}}]][t]==0.276281253576402052622079925762227839726037377161803163320316878162483075729117663799782079680408784210336126797460477138479369335828285`121.62355863066968,Derivative[1][\[DoubledGamma]][t]/4+Derivative[1][c[9]][t]==0,Derivative[1][\[DoubledGamma]][t]/5+Derivative[1][c[10]][t]==0,-(1/5) Derivative[1][\[DoubledGamma]][t]+E^10 Derivative[1][c[1]][t]==0,-Derivative[1][\[DoubledGamma]][t]+E^2 Derivative[1][c[5]][t]==0,-\[Lambda][{5,{9,5}}][t] (-Min[1,s[9,5][t]] Derivative[1][\[DoubledGamma]][t]+1/4 Min[s[5,9][t],s[9,5][t]] Derivative[1][\[DoubledGamma]][t]+s[5,9][t] Derivative[1][\[DoubledGamma]][t]+1/2 E^(2 Min[1,s[9,5][t]]) Derivative[1][c[5]][t]-1/2 E^(2 s[5,9][t]) Derivative[1][c[5]][t]+1/8 Derivative[1][c[9]][t]-1/8 E^(-8 Min[s[5,9][t],s[9,5][t]]) Derivative[1][c[9]][t])+\[Lambda][{9,{9,1}}][t] (1/4 Min[s[9,1][t],s[9,5][t]] Derivative[1][\[DoubledGamma]][t]+1/8 Derivative[1][c[9]][t]-1/8 E^(-8 Min[s[9,1][t],s[9,5][t]]) Derivative[1][c[9]][t])-\[Lambda][{5,{10,5}}][t] (-Min[1,s[9,5][t]] Derivative[1][\[DoubledGamma]][t]+1/5 Min[s[5,10][t],s[9,5][t]] Derivative[1][\[DoubledGamma]][t]+s[5,10][t] Derivative[1][\[DoubledGamma]][t]+1/2 E^(2 Min[1,s[9,5][t]]) Derivative[1][c[5]][t]-1/2 E^(2 s[5,10][t]) Derivative[1][c[5]][t]+1/10 Derivative[1][c[10]][t]-1/10 E^(-10 Min[s[5,10][t],s[9,5][t]]) Derivative[1][c[10]][t])-Derivative[1][k[5]][t]+Derivative[1][k[9]][t]-(3/2 Min[1,s[9,5][t]]+1/2 Min[1,s[9,5][t]]^2+33/32 Min[s[5,9][t],s[9,5][t]]-1/8 Min[s[5,9][t],s[9,5][t]]^2-Min[1,s[9,5][t]] \[DoubledGamma][t]+1/4 Min[s[5,9][t],s[9,5][t]] \[DoubledGamma][t]+1/2 E^(2 Min[1,s[9,5][t]]) c[5][t]-1/2 E^(2 s[5,9][t]) c[5][t]+1/8 c[9][t]-1/8 E^(-8 Min[s[5,9][t],s[9,5][t]]) c[9][t]-3/2 s[5,9][t]+\[DoubledGamma][t] s[5,9][t]-1/2 s[5,9][t]^2) Derivative[1][\[Lambda][{5,{9,5}}]][t]-(3/2 Min[1,s[9,5][t]]+1/2 Min[1,s[9,5][t]]^2+51/50 Min[s[5,10][t],s[9,5][t]]-1/10 Min[s[5,10][t],s[9,5][t]]^2-Min[1,s[9,5][t]] \[DoubledGamma][t]+1/5 Min[s[5,10][t],s[9,5][t]] \[DoubledGamma][t]+1/2 E^(2 Min[1,s[9,5][t]]) c[5][t]-1/2 E^(2 s[5,10][t]) c[5][t]+1/10 c[10][t]-1/10 E^(-10 Min[s[5,10][t],s[9,5][t]]) c[10][t]-3/2 s[5,10][t]+\[DoubledGamma][t] s[5,10][t]-1/2 s[5,10][t]^2) Derivative[1][\[Lambda][{5,{10,5}}]][t]+(33/32 Min[s[9,1][t],s[9,5][t]]-1/8 Min[s[9,1][t],s[9,5][t]]^2+1/4 Min[s[9,1][t],s[9,5][t]] \[DoubledGamma][t]+1/8 c[9][t]-1/8 E^(-8 Min[s[9,1][t],s[9,5][t]]) c[9][t]) Derivative[1][\[Lambda][{9,{9,1}}]][t]==0.276281253576402052622079925762227839726037377161803163320316878162483075729117663799782079680408784210336126797460477138479369335828285`121.62355863066968,-\[Lambda][{5,{9,5}}][t] (-Min[1,s[10,5][t]] Derivative[1][\[DoubledGamma]][t]+1/4 Min[s[5,9][t],s[10,5][t]] Derivative[1][\[DoubledGamma]][t]+s[5,9][t] Derivative[1][\[DoubledGamma]][t]+1/2 E^(2 Min[1,s[10,5][t]]) Derivative[1][c[5]][t]-1/2 E^(2 s[5,9][t]) Derivative[1][c[5]][t]+1/8 Derivative[1][c[9]][t]-1/8 E^(-8 Min[s[5,9][t],s[10,5][t]]) Derivative[1][c[9]][t])-\[Lambda][{5,{10,5}}][t] (-Min[1,s[10,5][t]] Derivative[1][\[DoubledGamma]][t]+1/5 Min[s[5,10][t],s[10,5][t]] Derivative[1][\[DoubledGamma]][t]+s[5,10][t] Derivative[1][\[DoubledGamma]][t]+1/2 E^(2 Min[1,s[10,5][t]]) Derivative[1][c[5]][t]-1/2 E^(2 s[5,10][t]) Derivative[1][c[5]][t]+1/10 Derivative[1][c[10]][t]-1/10 E^(-10 Min[s[5,10][t],s[10,5][t]]) Derivative[1][c[10]][t])+\[Lambda][{10,{10,1}}][t] (1/5 Min[s[10,1][t],s[10,5][t]] Derivative[1][\[DoubledGamma]][t]+1/10 Derivative[1][c[10]][t]-1/10 E^(-10 Min[s[10,1][t],s[10,5][t]]) Derivative[1][c[10]][t])-Derivative[1][k[5]][t]+Derivative[1][k[10]][t]-(3/2 Min[1,s[10,5][t]]+1/2 Min[1,s[10,5][t]]^2+33/32 Min[s[5,9][t],s[10,5][t]]-1/8 Min[s[5,9][t],s[10,5][t]]^2-Min[1,s[10,5][t]] \[DoubledGamma][t]+1/4 Min[s[5,9][t],s[10,5][t]] \[DoubledGamma][t]+1/2 E^(2 Min[1,s[10,5][t]]) c[5][t]-1/2 E^(2 s[5,9][t]) c[5][t]+1/8 c[9][t]-1/8 E^(-8 Min[s[5,9][t],s[10,5][t]]) c[9][t]-3/2 s[5,9][t]+\[DoubledGamma][t] s[5,9][t]-1/2 s[5,9][t]^2) Derivative[1][\[Lambda][{5,{9,5}}]][t]-(3/2 Min[1,s[10,5][t]]+1/2 Min[1,s[10,5][t]]^2+51/50 Min[s[5,10][t],s[10,5][t]]-1/10 Min[s[5,10][t],s[10,5][t]]^2-Min[1,s[10,5][t]] \[DoubledGamma][t]+1/5 Min[s[5,10][t],s[10,5][t]] \[DoubledGamma][t]+1/2 E^(2 Min[1,s[10,5][t]]) c[5][t]-1/2 E^(2 s[5,10][t]) c[5][t]+1/10 c[10][t]-1/10 E^(-10 Min[s[5,10][t],s[10,5][t]]) c[10][t]-3/2 s[5,10][t]+\[DoubledGamma][t] s[5,10][t]-1/2 s[5,10][t]^2) Derivative[1][\[Lambda][{5,{10,5}}]][t]+(51/50 Min[s[10,1][t],s[10,5][t]]-1/10 Min[s[10,1][t],s[10,5][t]]^2+1/5 Min[s[10,1][t],s[10,5][t]] \[DoubledGamma][t]+1/10 c[10][t]-1/10 E^(-10 Min[s[10,1][t],s[10,5][t]]) c[10][t]) Derivative[1][\[Lambda][{10,{10,1}}]][t]==0.276281253576402052622079925762227839726037377161803163320316878162483075729117663799782079680408784210336126797460477138479369335828285`121.62355863066968,\[Lambda][{5,{9,5}}][t] (Boole[s[5,9][t]<z[5,10][t]] (-Min[1,z[5,10][t]] Derivative[1][\[DoubledGamma]][t]+s[5,9][t] Derivative[1][\[DoubledGamma]][t]+1/2 E^(2 Min[1,z[5,10][t]]) Derivative[1][c[5]][t]-1/2 E^(2 s[5,9][t]) Derivative[1][c[5]][t])+Boole[0<z[5,10][t]] (1/4 Min[s[5,9][t],z[5,10][t]] Derivative[1][\[DoubledGamma]][t]+1/8 Derivative[1][c[9]][t]-1/8 E^(-8 Min[s[5,9][t],z[5,10][t]]) Derivative[1][c[9]][t]))+Boole[0<z[5,10][t]] \[Lambda][{5,{10,5}}][t] (1/5 Min[s[5,10][t],z[5,10][t]] Derivative[1][\[DoubledGamma]][t]+1/10 Derivative[1][c[10]][t]-1/10 E^(-10 Min[s[5,10][t],z[5,10][t]]) Derivative[1][c[10]][t])-Boole[0<z[5,10][t]] \[Lambda][{10,{10,1}}][t] (1/5 Min[s[10,1][t],z[5,10][t]] Derivative[1][\[DoubledGamma]][t]+1/10 Derivative[1][c[10]][t]-1/10 E^(-10 Min[s[10,1][t],z[5,10][t]]) Derivative[1][c[10]][t])+Derivative[1][k[5]][t]-Derivative[1][k[10]][t]+(Boole[0<z[5,10][t]] (33/32 Min[s[5,9][t],z[5,10][t]]-1/8 Min[s[5,9][t],z[5,10][t]]^2+1/4 Min[s[5,9][t],z[5,10][t]] \[DoubledGamma][t]+1/8 c[9][t]-1/8 E^(-8 Min[s[5,9][t],z[5,10][t]]) c[9][t])+Boole[s[5,9][t]<z[5,10][t]] (3/2 Min[1,z[5,10][t]]+1/2 Min[1,z[5,10][t]]^2-Min[1,z[5,10][t]] \[DoubledGamma][t]+1/2 E^(2 Min[1,z[5,10][t]]) c[5][t]-1/2 E^(2 s[5,9][t]) c[5][t]-3/2 s[5,9][t]+\[DoubledGamma][t] s[5,9][t]-1/2 s[5,9][t]^2)) Derivative[1][\[Lambda][{5,{9,5}}]][t]+Boole[0<z[5,10][t]] (51/50 Min[s[5,10][t],z[5,10][t]]-1/10 Min[s[5,10][t],z[5,10][t]]^2+1/5 Min[s[5,10][t],z[5,10][t]] \[DoubledGamma][t]+1/10 c[10][t]-1/10 E^(-10 Min[s[5,10][t],z[5,10][t]]) c[10][t]) Derivative[1][\[Lambda][{5,{10,5}}]][t]-Boole[0<z[5,10][t]] (51/50 Min[s[10,1][t],z[5,10][t]]-1/10 Min[s[10,1][t],z[5,10][t]]^2+1/5 Min[s[10,1][t],z[5,10][t]] \[DoubledGamma][t]+1/10 c[10][t]-1/10 E^(-10 Min[s[10,1][t],z[5,10][t]]) c[10][t]) Derivative[1][\[Lambda][{10,{10,1}}]][t]==0.276281253576402052622079925762227839726037377161803163320316878162483075729117663799782079680408784210336126797460477138479369335828285`121.62355863066968};

The differential equations to maintain the points s of intersection are:

sEqns = Flatten[Table[With[{i = i, j=j}, s[i,j]'[t]==(-(Derivative[1][\[DoubledGamma]][t]/u[i])+Derivative[1][\[DoubledGamma]][t]/u[j]-E^(-2 u[i] s[i,j][t]) Derivative[1][c[i]][t]+E^(-2 u[j] s[i,j][t]) Derivative[1][c[j]][t])/(-(1/u[i])+1/u[j]-2 E^(-2 u[i] s[i,j][t]) u[i] c[i][t]+2 E^(-2 u[j] s[i,j][t]) u[j] c[j][t])], {i, {1,5,9,10}}, {j, Complement[{1,5,9,10},{i}]}]]/.rates;

The equation to maintain z[5,9] as the point of intersection is:

zEqns ={Derivative[1][z[5,10]][t]==(Boole[s[5,9][t]<z[5,10][t]<=1] \[Lambda][{5,{9,5}}][t] Derivative[1][\[DoubledGamma]][t]-1/4 Boole[z[5,10][t]<=s[5,9][t]] \[Lambda][{5,{9,5}}][t] Derivative[1][\[DoubledGamma]][t]+Boole[s[5,10][t]<z[5,10][t]<=1] \[Lambda][{5,{10,5}}][t] Derivative[1][\[DoubledGamma]][t]-1/5 Boole[z[5,10][t]<=s[5,10][t]] \[Lambda][{5,{10,5}}][t] Derivative[1][\[DoubledGamma]][t]-1/5 Boole[s[10,1][t]<z[5,10][t]<=1] \[Lambda][{10,{10,1}}][t] Derivative[1][\[DoubledGamma]][t]+1/5 Boole[z[5,10][t]<=s[10,1][t]] \[Lambda][{10,{10,1}}][t] Derivative[1][\[DoubledGamma]][t]+E^(10 z[5,10][t]) Boole[s[10,1][t]<z[5,10][t]<=1] \[Lambda][{10,{10,1}}][t] Derivative[1][c[1]][t]-E^(2 z[5,10][t]) Boole[s[5,9][t]<z[5,10][t]<=1] \[Lambda][{5,{9,5}}][t] Derivative[1][c[5]][t]-E^(2 z[5,10][t]) Boole[s[5,10][t]<z[5,10][t]<=1] \[Lambda][{5,{10,5}}][t] Derivative[1][c[5]][t]-E^(-8 z[5,10][t]) Boole[z[5,10][t]<=s[5,9][t]] \[Lambda][{5,{9,5}}][t] Derivative[1][c[9]][t]-E^(-10 z[5,10][t]) Boole[z[5,10][t]<=s[5,10][t]] \[Lambda][{5,{10,5}}][t] Derivative[1][c[10]][t]+E^(-10 z[5,10][t]) Boole[z[5,10][t]<=s[10,1][t]] \[Lambda][{10,{10,1}}][t] Derivative[1][c[10]][t]-(Boole[z[5,10][t]<=s[5,9][t]] (33/32+\[DoubledGamma][t]/4+E^(-8 z[5,10][t]) c[9][t]-1/4 z[5,10][t])+Boole[s[5,9][t]<z[5,10][t]<=1] (3/2-\[DoubledGamma][t]+E^(2 z[5,10][t]) c[5][t]+z[5,10][t])) Derivative[1][\[Lambda][{5,{9,5}}]][t]-(Boole[z[5,10][t]<=s[5,10][t]] (51/50+\[DoubledGamma][t]/5+E^(-10 z[5,10][t]) c[10][t]-1/5 z[5,10][t])+Boole[s[5,10][t]<z[5,10][t]<=1] (3/2-\[DoubledGamma][t]+E^(2 z[5,10][t]) c[5][t]+z[5,10][t])) Derivative[1][\[Lambda][{5,{10,5}}]][t]+(Boole[z[5,10][t]<=s[10,1][t]] (51/50+\[DoubledGamma][t]/5+E^(-10 z[5,10][t]) c[10][t]-1/5 z[5,10][t])+Boole[s[10,1][t]<z[5,10][t]<=1] (51/50-\[DoubledGamma][t]/5+E^(10 z[5,10][t]) c[1][t]+1/5 z[5,10][t])) Derivative[1][\[Lambda][{10,{10,1}}]][t])/(Boole[s[5,9][t]<z[5,10][t]<=1] \[Lambda][{5,{9,5}}][t]-1/4 Boole[z[5,10][t]<=s[5,9][t]] \[Lambda][{5,{9,5}}][t]+2 E^(2 z[5,10][t]) Boole[s[5,9][t]<z[5,10][t]<=1] c[5][t] \[Lambda][{5,{9,5}}][t]-8 E^(-8 z[5,10][t]) Boole[z[5,10][t]<=s[5,9][t]] c[9][t] \[Lambda][{5,{9,5}}][t]+Boole[s[5,10][t]<z[5,10][t]<=1] \[Lambda][{5,{10,5}}][t]-1/5 Boole[z[5,10][t]<=s[5,10][t]] \[Lambda][{5,{10,5}}][t]+2 E^(2 z[5,10][t]) Boole[s[5,10][t]<z[5,10][t]<=1] c[5][t] \[Lambda][{5,{10,5}}][t]-10 E^(-10 z[5,10][t]) Boole[z[5,10][t]<=s[5,10][t]] c[10][t] \[Lambda][{5,{10,5}}][t]-1/5 Boole[s[10,1][t]<z[5,10][t]<=1] \[Lambda][{10,{10,1}}][t]+1/5 Boole[z[5,10][t]<=s[10,1][t]] \[Lambda][{10,{10,1}}][t]-10 E^(10 z[5,10][t]) Boole[s[10,1][t]<z[5,10][t]<=1] c[1][t] \[Lambda][{10,{10,1}}][t]+10 E^(-10 z[5,10][t]) Boole[z[5,10][t]<=s[10,1][t]] c[10][t] \[Lambda][{10,{10,1}}][t])};

Without "StiffnessSwitching" NDSolve is only able to reduce t to 0.991 in 500 steps. With "StiffnessSwitching" NDSolve successfully reduces t to the point at which [Lambda][{5, {9,5}}][t] reaches 0 and the system of equations changes in fewer than 500 steps.

With[{opts=SystemOptions[]}, Off[SystemOptions::noset]; SetSystemOptions["NDSolveOptions"->"DefaultSolveTimeConstraint"->5000.`]; sol2 = First[NDSolve[Rationalize[Join[initvals, mainEqns, sEqns, zEqns, {WhenEvent[\[Lambda][{5,{9,5}}][t]<0,"StopIntegration"]}],0], Join[mainVars, sVars, zVars], {t, 0, 1}, MaxSteps ->500, WorkingPrecision->60, Method->"StiffnessSwitching"]]; SetSystemOptions[opts];]

The dialog at the top shows how the solution changes with t.(Edit: I was unable to get iframe embedded in the question: Here's a link instead: How solution changes with t) Reduce the "max x" to 0.2 and observe how ```[DoubledGamma]`` reduces and f[x,5,t] moves between f[x,9,t] and f[x,10,t] on the domain [0, s[5,9][t]] as you reduce t.

Variable Names Matter Simply changing \[DoubledGamma]->obj and \[Lambda]->w leads to Power::infy error. This is the same behavior exhibited with the original names when solving in Mathematica Versions 14.3.0.0, 14.2.1.0 and 14.2.0.0.

With[{opts=SystemOptions[]}, Off[SystemOptions::noset]; SetSystemOptions["NDSolveOptions"->"DefaultSolveTimeConstraint"->5000.`]; sol2 = First[NDSolve[Rationalize[Join[initvals, mainEqns, sEqns, zEqns, {WhenEvent[\[Lambda][{5,{9,5}}][t]<0,"StopIntegration"]}],0]/.\[DoubledGamma]->obj/.\[Lambda]->w, Join[mainVars, sVars, zVars]/.\[DoubledGamma]->obj/.\[Lambda]->w, {t, 0, 1}, MaxSteps ->500, WorkingPrecision->60, Method->"StiffnessSwitching"]]; SetSystemOptions[opts];]

Error Messages:

Power::infy: Infinite expression 1/0.*10^-54 encountered.

Power::infy: Infinite expression 1/0.*10^-63 encountered.

Infinity::indet: Indeterminate expression 0.*10^-124 E^14 ComplexInfinity ComplexInfinity encountered.

Power::infy: Infinite expression 1/0.*10^-54 encountered.

General::stop: Further output of Power::infy will be suppressed during this calculation.

Infinity::indet: Indeterminate expression 0.*10^-68 E^22 ComplexInfinity ComplexInfinity encountered.

Infinity::indet: Indeterminate expression 0.*10^-124 E^4 ComplexInfinity ComplexInfinity encountered.

General::stop: Further output of Infinity::indet will be suppressed during this calculation.

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 1.`60..

Improve this question
$\endgroup$

1 Answer 1

Reset to default
-1
$\begingroup$

My initial answer appears in retrospect to be incorrect and was kindly deleted. This answer resolves the issue, but is considerably slower. The issue arises when NDSolve`EquationProcessing transforms the system of equations into explicit form. Indeed, the problem occurs even, recognizing that once we solve mainEqns for the mainVars, applying that solution to sEqns and zEqns expresses them in explicit form.

The following approaches, based on the knowledge that the last [DoubledGamma] variable and the last mainEqn are special, failed because they took too long (hours in some cases) even for this small example:

  1. Eliminate the last [DoubledGamma] from mainEqns using Eliminate. Eliminate takes a very long time.
  2. Construct explicit definitions sequentially by a. Solving all but the last mainEqn for all but the last last \[DoubledGamma] variable and then b. Transforming the last mainEqn with that solution and solving the resulting equation for the \[DoubledGamma] variable. That produced explicit expressions in reasonable time, but NDSolve`EquationProcessing takes a very long time to process the definitions even with "SystemSimplification" ->False.
  3. Solve with Assumptions, e.g., that z[5,10][t]==0 and z[5,10][t]<s[5,9][t]. With those Assumptions, Solve takes a very long time.

What finally worked is to construct explicit solutions using LinearSolve. That takes, surprisingly, about twice as long as Solve to construct the explicit expressions, but produces the correct ones. To ensure robustness, I added ZeroTest ->(N[#/.(initvals/.Equal->Rule/.exp__[1]:>exp[t]), 60 (*The selected WorkingPrecision*)]==0&) to avoid division by an expression that is 0 under the initial values. Adding that test had little impact on the timing.

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.