
Mathematica: Initial value trouble when using NDSolve for differentialalgebraic system
Posted:
Oct 19, 2011 1:01 PM


Hello to all!
I would greatly appreciate your suggestions about resolving a problem with initial conditions for differentialalgebraic 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}}

