# 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] == 0, r[0] == 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] == 0, y[0] == 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

1
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

1
@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

1
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
show 8 more comments

## 2 Answers

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.

OLD ANSWER

(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[0] \[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.♦ 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] == 0,
r[0] == 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
lang-mma
{"extraction_info": {"found_math": true, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 28, "mathjax_display_tex": 0, "mathjax_asciimath": 0, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0}, "config": {"markdown_headings": true, "markdown_code": true, "boilerplate_config": {"ratio_threshold": 0.18, "absolute_threshold": 10, "end_threshold": 15, "enable": true}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "github"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "english_confidence": 0.9112802147865295, "perplexity_flag": "middle"}