Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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"];
share|improve this question
    
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

up vote 1 down vote accepted

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
 ]

Mathematica graphics

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]"}
 ]

Mathematica graphics

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.

share|improve this answer
    
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
add comment

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"]

enter image description here

share|improve this answer
    
Well, no, the problem is that I don't know what x should be. –  Andrew May 19 '12 at 14:18
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.