Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

Views expressed in these public forums are not endorsed by Drexel University or The Math Forum.


Math Forum » Discussions » sci.math.* » sci.math.symbolic.independent

Topic: Mathematica: Initial value trouble when using NDSolve for
differential-algebraic system

Replies: 1   Last Post: Oct 26, 2011 2:25 AM

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
Alexander Brown

Posts: 2
Registered: 10/19/11
Mathematica: Initial value trouble when using NDSolve for
differential-algebraic system

Posted: Oct 19, 2011 1:01 PM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

Hello to all!

I would greatly appreciate your suggestions about resolving a
problem with initial conditions for differential-algebraic system in
Mathematica.
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}}



Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.