Prev: V. 6.0: PNG export problems
Next: Numerical solution of differential equation
From: Zachar István on 22 May 2010 00:40
Dear Group,

I have an ODE set, which I send in to NDSolve. It solve more-or-less
correctly, but indicates a stiffness somewhere. As a matter of fact, I
cannot rely on the InterpolatingFunction output, so I'm saving the
actual datapoints via EvaluationMonitor (or StepMonitor). But the
monitor stops when the stiffness is found, while the
InterpolatingFunction goes on. Why does this happen? How can I collect
ALL the datapoints from the integration, to produce a ListLinePlot
identical to the Plot output?
Any help would be appreciated.

eq = {
Derivative[1][A1][t] == -0.2*A1[t]^2 + 0.1*A2[t] + 20*A5[t] -
2*A1[t]*X[t],
Derivative[1][A2][t] == -100.1*A2[t] + 2*A1[t]*X[t] +
0.1*A3[t]*Y[t],
Derivative[1][A3][t] ==
100*A2[t] - 100*A3[t] + 0.1*A4[t]*U[t] - 0.1*A3[t]*Y[t],
Derivative[1][A4][t] ==
100*A3[t] - 100*A4[t] + 0.1*A5[t]*Tp[t] - 0.1*A4[t]*U[t],
Derivative[1][A5][t] ==
0.1*A1[t]^2 + 100*A4[t] - 10*A5[t] - 0.1*A5[t]*Tp[t],
Derivative[1][d][t] == X[t] + Y[t] + Z1[t] + Z2[t],
Derivative[1][pV0][t] ==
0.1*pV1[t]*R[t] - pV0[t]*V[t] + 2*pV5[t]*V[t],
Derivative[1][pV1][t] == -0.1*pV1[t]*R[t] + pV0[t]*V[t] -
pV1[t]*V[t], Derivative[1][pV2][t] == pV1[t]*V[t] - pV2[t]*V[t],
Derivative[1][pV3][t] == pV2[t]*V[t] - pV3[t]*V[t],
Derivative[1][pV4][t] == pV3[t]*V[t] - pV4[t]*V[t],
Derivative[1][pV5][t] == pV4[t]*V[t] - pV5[t]*V[t],
Derivative[1][pW0][t] ==
0.1*pW1[t]*R[t] - pW0[t]*W[t] + 20*pW5[t]*W[t],
Derivative[1][pW1][t] == -0.1*pW1[t]*R[t] + pW0[t]*W[t] -
10*pW1[t]*W[t],
Derivative[1][pW2][t] == 10*pW1[t]*W[t] - 10*pW2[t]*W[t],
Derivative[1][pW3][t] == 10*pW2[t]*W[t] - 10*pW3[t]*W[t],
Derivative[1][pW4][t] == 10*pW3[t]*W[t] - 10*pW4[t]*W[t],
Derivative[1][pW5][t] == 10*pW4[t]*W[t] - 10*pW5[t]*W[t],
Derivative[1][R][t] == -0.1*pV1[t]*R[t] - 0.1*pW1[t]*R[t] +
0.1*T[t] - 10*R[t]*Ts[t] + pV0[t]*V[t] + pV1[t]*V[t] +
pV2[t]*V[t] + pV3[t]*V[t] + pV4[t]*V[t] + pV5[t]*V[t] +
pW0[t]*W[t] + 10*pW1[t]*W[t] + 10*pW2[t]*W[t] + 10*pW3[t]*W[t] +
10*pW4[t]*W[t] + 10*pW5[t]*W[t],
Derivative[1][T][t] == -10.1*T[t] + 10*R[t]*Ts[t],
Derivative[1][Tm][t] == 10*T[t],
Derivative[1][Tp][t] == 100*A4[t] - 10*Tp[t] - 0.1*A5[t]*Tp[t],
Derivative[1][Ts][t] == 0.1*T[t] + 10*Tp[t] - 10*R[t]*Ts[t],
Derivative[1][U][t] ==
100*A3[t] - 0.1*A4[t]*U[t] + 0.1*V[t] + 0.1*W[t] -
1.*U[t]*Z1[t] - 1.*U[t]*Z2[t],
Derivative[1][V][t] ==
0.1*pV1[t]*R[t] - 0.1*V[t] - pV0[t]*V[t] - pV1[t]*V[t] -
pV2[t]*V[t] - pV3[t]*V[t] - pV4[t]*V[t] - pV5[t]*V[t] +
1.*U[t]*Z1[t],
Derivative[1][W][t] ==
0.1*pW1[t]*R[t] - 0.1*W[t] - pW0[t]*W[t] - 10*pW1[t]*W[t] -
10*pW2[t]*W[t] - 10*pW3[t]*W[t] - 10*pW4[t]*W[t] -
10*pW5[t]*W[t] + 1.*U[t]*Z2[t],
Derivative[1][X][t] == 20 + 0.1*A2[t] - X[t] - 2*A1[t]*X[t],
Derivative[1][Y][t] == 100*A2[t] - Y[t] - 0.1*A3[t]*Y[t],
Derivative[1][Z1][t] == 10 + 0.1*V[t] - Z1[t] - 1.*U[t]*Z1[t],
Derivative[1][Z2][t] == 10 + 0.1*W[t] - Z2[t] - 1.*U[t]*Z2[t],
A1[0] == 5, A2[0] == 0, A3[0] == 0, A4[0] == 0, A5[0] == 0,
d[0] == 0, pV0[0] == 1, pV1[0] == 0, pV2[0] == 0, pV3[0] == 0,
pV4[0] == 0, pV5[0] == 0, pW0[0] == 0, pW1[0] == 0, pW2[0] == 0,
pW3[0] == 0, pW4[0] == 0, pW5[0] == 0, R[0] == 0, T[0] == 0,
Tm[0] == 50, Tp[0] == 0, Ts[0] == 0, U[0] == 0, V[0] == 0,
W[0] == 0, X[0] == 0, Y[0] == 0, Z1[0] == 0, Z2[0] == 0};

all = {A1, A2, A3, A4, A5, d, pV0, pV1, pV2, pV3, pV4, pV5, pW0, pW1,
pW2, pW3, pW4, pW5, R, T, Tm, Tp, Ts, U, V, W, X, Y, Z1, Z2};

time = 20;

data = Reap[NDSolve[eq, all, {t, 0, time},
EvaluationMonitor :> Sow[Evaluate[{{t, pV0[t]}, {t, pW0[t]}}]],
Method -> "ExplicitRungeKutta"
]]

Transpose(a)First@Last(a)data;

Plot[Evaluate[{pV0[t], pW0[t]} /. First(a)data], {t, 0, time}, ImageSize
-> 500, Frame -> True, PlotRange -> {{0, time}, {0, 7}},
PlotRangePadding -> .5]

ListLinePlot[Transpose(a)First@Last(a)data, ImageSize -> 500, Frame ->
True, PlotRange -> {{0, time}, {0, 7}}, PlotRangePadding -> .5]

 | 
Pages: 1
Prev: V. 6.0: PNG export problems
Next: Numerical solution of differential equation