MathGroup Archive 2011

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Initial value problem when using NDSolve for differential-algebraic system

  • To: mathgroup at smc.vnet.net
  • Subject: [mg122228] Re: Initial value problem when using NDSolve for differential-algebraic system
  • From: Oliver Ruebenkoenig <ruebenko at wolfram.com>
  • Date: Fri, 21 Oct 2011 06:25:32 -0400 (EDT)
  • Delivered-to: l-mathgroup@mail-archive0.wolfram.com

Alexander,

you could try to use Rationalize on you input and derive the necessary 
precision from that or have NDSolve do it then.

Hope this gets you started.

Oliver

On Thu, 20 Oct 2011, Alexander Brown wrote:

> Hello to all!
>
>  I would greatly appreciate your suggestions about resolving a
> problem with initial conditions for differential-algebraic system.
>  The system is BelModelTanhSimpleAlien (see below). Since NDSolve is
> known to have trouble with finding initial conditions for its
> variables, I used Solve to find a solution for algebraic part
> (solFullEquilibrium) of my equations at t=0 and then used the values
> obtained (jcPKA, ycAMP, yCPKI) as initial values for NDSolve. Rather
> annoyingly, I still got the icfail error from NDSolve. Then,
> naturally, I've checked whether (jcPKA, ycAMP, yCPKI) actually solves
> solFullEquilibrium. Substitution gave {False} for all equations, with
> the left equations and difference between left and right sides being
> at the order of 10^-3 and 10^-17, respectively.
>
>  I think that icfail error happens because NDSolve uses too much
> precision when checking initial values, but I don't know how to lower
> it (actually, can't lower the precision for solFullEquilibrium,
> either). I would appreciate your suggestions to deal with this
> problem.
>
>  You could find the code for searching initial values (1), trying to
> launch NDSolve (2), and checking the obtained solution (3) below.
>
> Thank you in advance,
> Alexander Brown
>
> 1.Solve input (search initial values):
>
> eqFullEquilibrium
> solFullEquilibrium =
> Solve[Flatten[{eqFullEquilibrium}], {jcPKA, ycAMP, yCPKI}, Reals]
>
> Solve output:
>
> {0.0015 == (1 +
>     228.571 jcPKA (1 + (3280. (1 + 18280./(1 + Tanh[ycAMP])))/(
>        1 + Tanh[ycAMP]))) (jcPKA + 0.00015 (1 + Tanh[yCPKI])),
> 1/1000000 == (1 + Tanh[ycAMP])/2000000 +
>   2 (1 + 228.571 jcPKA (1 + 1640./(1 + Tanh[ycAMP]))) (jcPKA +
>      0.00015 (1 + Tanh[yCPKI])),
> 0.00015 (1 + Tanh[yCPKI]) ==
>  10000. jcPKA (0.0003 - 0.00015 (1 + Tanh[yCPKI])),
> jcPKA >= 0, (1 + Tanh[ycAMP])/2000000 >= 0,
> 0.00015 (1 + Tanh[yCPKI]) >= 0}
>
> {{jcPKA -> 8.68594*10^-8, ycAMP -> -0.516519, yCPKI -> -3.52432}}
>
> NDSolve input (try to launch NDSolve):
>
> BelModelTanhSimpleAlien
> BelModelTanhSimpleAlienVar
> lsol = NDSolve[{BelModelTanhSimpleAlien, xcAMP[0] ==
> -0.5165187286029376`, xCPKI[0] == -3.524317593645514`},
> BelModelTanhSimpleAlienVar,
> {t, 0, 20},
> MaxSteps -> Infinity, MaxStepFraction -> 0.0001, AccuracyGoal -> 2,
> PrecisionGoal -> 2]
>
> 2.NDSolve output:
>
> {Derivative[1][cPKA][
>   t] == -0.03 cPKA[t] (1 + 1/2 (-1 - Tanh[xCPKI[t]])) +
>   1.5*10^-6 (1 + Tanh[xCPKI[t]]),
> 0.0015 == (1 +
>     228.571 cPKA[
>       t] (1 + (
>        0.00328 (1 + 0.01828/(cAMP0[t] (1 + Tanh[xcAMP[t]]))))/(
>        cAMP0[t] (1 + Tanh[xcAMP[t]])))) (cPKA[t] +
>     0.00015 (1 + Tanh[xCPKI[t]])),
> cAMP0[t] ==
>  1/2 cAMP0[t] (1 + Tanh[xcAMP[t]]) +
>   2 (1 + 228.571 cPKA[
>        t] (1 + 0.00164/(cAMP0[t] (1 + Tanh[xcAMP[t]])))) (cPKA[t] +
>      0.00015 (1 + Tanh[xCPKI[t]])),
> Derivative[1][PP][t] ==
>  0.1 (1/1000 - PP[t]) - (10 cPKA[t] PP[t])/(3/200 + PP[t]),
> Derivative[1][PhKa][t] == (30 cPKA[t] (1/200 - PhKa[t]))/(
>   1/100 - PhKa[t]) - (0.5 PhKa[t] PP[t])/(1/250 + PhKa[t]),
> Derivative[1][GPa][t] == (30 (7/100 - GPa[t]) PhKa[t])/(
>   3/40 - GPa[t]) - (12 GPa[t] PP[t])/(1/250 + GPa[t]),
> cPKA[0] == 8.68594*10^-8, PP[0] == 0.000999457,
> PhKa[0] == 0.0000104451, GPa[0] == 0.0000999686,
> cAMP0[t] ==
>  1/1000000 + 0.001 2^(-(1/4) (-730 + t)^2) +
>   0.002 2^(-(1/4) (-630 + t)^2) + 0.003 2^(-(1/4) (-530 + t)^2) +
>   0.004 2^(-(1/4) (-430 + t)^2) + 0.005 2^(-(1/4) (-330 + t)^2) +
>   0.02 2^(-(1/4) (-230 + t)^2) + 0.4 2^(-(1/4) (-130 + t)^2) +
>   5 2^(-(1/4) (-30 + t)^2)}
>
> {cAMP0, xcAMP, cPKA, xCPKI, PP, PhKa, GPa}
>
> NDSolve::icfail: Unable to find initial conditions that satisfy the
> residual function within specified tolerances. Try giving initial
> conditions for both values and derivatives of the functions. >>
>
> {}
>
>
> 3.Check solution from Solve input:
>
> eqFullEquilibrium /. solFullEquilibrium
> eqFullEquilibrium /. Equal -> Subtract /. solFullEquilibrium
> (*The following is an attempt to lower precision. Didn't work.*)
> Block[{$MinPrecision = 2, $MaxPrecision = 2}, eqFullEquilibrium /.
> solFullEquilibrium]
> Block[{$MinPrecision = 2, $MaxPrecision = 2}, eqFullEquilibrium /.
> Equal -> Subtract /. solFullEquilibrium]
>
> Check solution from Solve output:
>
> {{False, False, False, True, True, True}}
> {{-3.23092*10^-17, -5.39984*10^-20, 7.30566*10^-21, True, True, True}}
> {{False, False, False, True, True, True}}
> {{-3.23092*10^-17, -5.39984*10^-20, 7.30566*10^-21, True, True, True}}
>
>

-- 



  • Prev by Date: Re: Problem with "Which"
  • Next by Date: A couple of questions
  • Previous by thread: Initial value problem when using NDSolve for differential-algebraic system
  • Next by thread: Single precision floating point numbers