Message Boards Message Boards

GROUPS:

NDSolve error when using an oscillatory forcing of ODE system.

Posted 3 months ago
407 Views
|
3 Replies
|
0 Total Likes
|

I am experiencing problems when using NDSolve to find the time-dependent solution to a system of 37 coupled ODEs. ODEs contain terms which are forced through oscillatory sine functions.

This takes place in a numerical model where variables for which I want to find a time-dependent solution are the concentration of chemical species in an ocean chemistry model, while model forcing is the oscillatory evolution of the hydrologcal cycle (evaporation and river runoff water fluxes).

The code is very long and I will not copy it here, hoping that this text will be enough for some of you to guide me in a promising direction.

(1) My base setting of the NDSolve function was:

res = NDSolve[Join[ODEsys, IC], var, {t, 0, tmax}, AccuracyGoal -> 14,
   PrecisionGoal -> 14]

=> this returned the following error:

NDSolve::ndcf: Repeated convergence test failure at t == 2.4122636505104266`*^7; unable to continue.

(2) I was able to obtain a solution for simulated times tmax smaller than 10 years using the following setting:

res = NDSolve[Join[ODEsys, IC], var, {t, 0, tmax}, AccuracyGoal -> 14,
   PrecisionGoal -> 14, 
  Method -> {"EquationSimplification" -> "Residual"}, 
  MaxStepSize -> tmax/150000, MaxSteps -> Infinity]

(3) However, for longer simulation times (tmax > 10 years), I get the following error:

NDSolve::ndsz: At t == 2.9634156107275105`*^9, step size is effectively zero; singularity or stiff system suspected.

(4) Switching to the following setting which in other occasions solved stiff systems:

res = NDSolve[Join[ODEsys, IC], var, {t, 0, tmax}, AccuracyGoal -> 16,
   PrecisionGoal -> 16, WorkingPrecision -> 18, 
  Method -> {"StiffnessSwitching"}, MaxStepSize -> tmax/150000, 
  MaxSteps -> Infinity]

produces this error message:

NDSolve::precw: The precision of the differential equation ({{(TemporaryVariable$1623^\[Prime])[t]==0,<<49>>,<<24>>},{},{},{},{}}) is less than WorkingPrecision (18.`). Tensors....have incompatible shapes.

I am stuck since my knowledge of the numerous settings of the NDSolve function is very limited.

Can anyone give me a tip or direct me towards documentation where i can learn more about how to solve these kinds of issues with NDSolve ?

Thank you for having read to this point!

Vanni

3 Replies
Posted 3 months ago

Crossposted here.

Giovanni,

As to your last error, you can't set precision to 16 or greater unless you have no machine precision numbers in your equations. You would need to make sure your coefficients are higher precision. This tutorial will explain the issue. It's hard to say more without an example.

Regards,

Neil

Dear Neal, thank you for your comment. I have now set the precision of all of my coefficients to 20 with SetPrecision[x,20], where x is a given coefficient. This indeed overcame problem nr. 4 in my message above. However, I am still left with an unsolvable stiff system, even if I impose WorkingPrecision->19. I will try increasing further the precision of my coefficients and setting a higher working precision. Also, I can try to modify my equation system, I suspect there might be steep chages in some derivatives that I can try to avoid. Cheers, Vanni

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract