NDSolve and results of the previous computations

I want to solve following system of ODEs:

$\Bigg\{ \begin{array}{} \frac{\partial C}{\partial t}=\frac{2W_b}{Rsin(\theta)}(1-\frac{C}{\gamma})+\frac{\nu}{\pi R^2}\frac{\partial C}{\partial Z} \\ \frac{\partial R}{\partial t}=\frac{W_b}{\gamma sin(\theta)} \end{array}$

Where $W_b=k_bC_h(1-\frac{C}{C_h})^{\frac{4}{3}}$ and $\nu, k_b, C_h, h, \gamma, \theta$ are constants.

It is known that:

$\frac{\partial C}{\partial Z}=\frac{C_i-C_{i-1}}{Z_i-Z_{i-1}}$, where $Z_i-Z_{i-1}=h$ and $h$ is a constant.

So, the first equation become this: $\frac{\partial C}{\partial t}=\frac{2W_b}{Rsin(\theta)}(1-\frac{C}{\gamma})+\frac{\nu}{\pi R^2}\frac{C_i-C_{i-1}}{h}$

My question is how can I get the result of previous computation, namely $C_{i-1}$?

Here is the mathematica code I wrote:

Wb[C_] := kb*Ch*(1 - C/Ch)^(4/3);
system := {
c'[t] == (2*Wb[c[t]])/(r[t]*Sin[theta])*(1 - C/gamma) + v/(Pi*r[t]^2)*(c[t] - ?? )/h,
r'[t] == Wb[c[t]]/(gamma*Sin[theta]),
c == 0, r == 0.15};
solution = First@NDSolve[system, {c, r}, {t, 0, 25900000}, Method -> "BDF"];

UPDATE: Simplified version of the problem.

This is a model based on a plug flow reactor.
$\Bigg\{ \begin{array}{} \dot{x} = \frac{c_1 x f(x)}{y} + \frac{\partial x}{\partial z} \\ \dot{y} = c_2 f(x) \end{array}$

Where c1 and c2 are constants. This is a model of a physical process and it was shown that z dimension is quantified by a chunks of a constant size h and x is monotonously increasing, so $\frac{\partial x}{\partial z}=\frac{\Delta x}{\Delta z}=\frac{x_i - x_{i-1}}{h}$

f[x_] := ...;
system := {
x'[t] == c1*x*f[x[t]]/y[t] + (x[t] - ?)/h,
y'[t] == c2*f[x[t]],
x == 0, y == 0.15};
solution = First@NDSolve[system, {x, y}, {t, 0, 25900000}, Method -> "BDF"];
• First of all, you can already solve for R(t) by hand. Second, the remaining equation probably should have C replaced by C_i. Is that what you mean? You will then also need the initial conditions for all C_i, I would guess. With that, you would have a coupled system of first-order equations in time for the C_i. – Jens May 19 '12 at 4:50
• I've updated the question. Wb is a function dependent on C, I'm not sure if it possible to solve it by hand. – Andrew May 19 '12 at 14:20
• So, initially $x=x(y,z)$ but you would like to replace this 2d function by a finite number of 1d functions, ie, you'd like to discretise the $z$ direction by hand. Right? – acl May 19 '12 at 17:07
• @Andrew then you need to solve a set of differential equations for $x_i(y,t)$ (ie, a set of variables, not a single var) and $y(x_1,\ldots,x_N,t)$, whereas you're trying to set the problem up for a single $x$. – acl May 19 '12 at 17:30
• Look up Delayed Differential Equations in the documentation (howto/SolveDelayDifferentialEquations). Your ?? will be x[t-h]. Also figure out the correct initial condition because you will need to specify it for a range of values. – Daniel Lichtblau May 20 '12 at 16:10

It seems (from the comments) that what you want to do is this: initially x=x(y,z) but you would like to replace this function of 2 vars by a finite number of functions of 1 var, ie, you'd like to discretize the $z$ direction by hand. This means that you need to solve a set of differential equations for the $x_i(y,t)$ (a set of functions, not a single fun) and $y(x_1,\ldots,x_N,t)$.

You're trying to set the problem up for a single $x$, and that is the problem.

(I think this is neat so I'll leave it here for now)

I very likely misunderstood you. If, however, I understood correctly the question, it is something like: Suppose I have $y'(t)=f(y)$ and ask mathematica to solve it numerically. How do I inspect the values of $t$ used?

The answer is to rig the ODE so that you can peek at the values mathematica evaluates, as follows:

ClearAll[f];
upT = 2*Pi;
f[y_?NumericQ, t_?NumericQ] := (Sow[{t, y}]; -Sin@y)

Here I define $f(y)=-\sin(y)$, but in a way that allows me to watch which values are passed to $f$ by NDSolve (the Sow bit). Then:

points = (sol = NDSolve[
{
y'[t] == f[y[t], t],
y \[Equal] 1
},
y,
{t, 0, upT}
]) // Reap // Last // Last;

solves the ODE, collecting the values of $y$ and $t$ passed to $f$ by NDSolve. Then plot the solution, along with the points at which $f$ was evaluated:

Plot[
y[t] /. sol,
{t, 0, upT},
Epilog :> {Red, PointSize[.015], Point[points]},
ImageSize -> 640
] Not sure if I am answering the right question though.

It's also interesting to look at the size of the steps taken:

ListPlot[
Differences[points[[All, 1]]],
ImageSize -> 640,
BaseStyle -> FontSize -> 20,
AxesLabel -> {"i", "t[i]-t[i-1]"}
] so they are not all equal (not all $h$), sometimes the solver goes backwards etc. Of course this will depend on the solver you use.

• You know, you could have used EvaluationMonitor :> Sow[{t, y[t]}] as an NDSolve[] option setting, don't you? – J. M. will be back soon May 19 '12 at 18:32
• @J.M. yes, but I am not the most elegant of programmers :) – acl May 19 '12 at 18:35

Your syntax has me a bit confused. SetDelayed (:=) looks off to me. Typically one uses this to define functions. I think you want Equal (==) where you have the system of equations and Set (=) where you assign the value of the expression to "solution".

Does this get you any closer (I replaced your "??" with "x")?

solution = First@NDSolve[{
c'[t] == (2*Wb[c[t]])/(r[t]*Sin[theta])*(1 - C/gamma) +
v/(Pi*r[t]^2)*(c[t] - x)/h,
r'[t] == Wb[c[t]]/(gamma*Sin[theta]),
c == 0,
r == 0.15
}, {c, r}, {t, 0, 25900000}, Method -> "BDF"] • Well, no, the problem is that I don't know what x should be. – Andrew May 19 '12 at 14:18