A system of ODEs (Point Kinetics)

Question

I seem to have reached a dead end while attempting to solve a system of ODEs that is supposed to model a point-kinetic reactor.

The situation is to model a ramp insertion of reactivity. The constants $a$, $\beta$, $\lambda$, $\Lambda$ are all parameters with given numerical values.

Here is my work thus far:

$$\begin{split} \frac{\mathrm{d}n(t)}{\mathrm{d}t} &= \frac{\rho(t) - \beta}{\Lambda} n(t) - \lambda c(t) \\ \frac{\mathrm{d}c(t)}{\mathrm{d}t} &= \frac{\beta}{\Lambda} n(t) + \lambda c(t) \\ \rho(t) &= at \end{split} $$

Then we take the Laplace Transform:

$$\begin{split} -n(0) + sN(s) &= - \frac{a}{\Lambda} \frac{\mathrm{d}N(s)}{\mathrm{d}s} - \frac{\beta}{\Lambda} N(s) + \lambda C(s) \\ -c(0) + sC(s) &= \frac{\beta}{\Lambda} - \lambda C(s) \\ \end{split}$$

Because the second equation is now no longer an ODE, we can just solve for $C(s)$ and plug it into the first:

$$ C(s) = \frac{\frac{\beta}{\Lambda} N(s) + c(0)}{s - \lambda} $$

$$ -n(0) + sN(s) = - \frac{a}{\Lambda} \frac{\mathrm{d}N(s)}{\mathrm{d}s} - \frac{\beta}{\Lambda} N(s) + \frac{\lambda \beta}{\Lambda(s + \lambda)} N(s) + \frac{\lambda c(0)}{s + \lambda} $$

Then place it into standard form:

$$ \frac{\mathrm{d}N(s)}{\mathrm{d}s} + \frac{\Lambda}{a}\left(s + \frac{\beta}{\Lambda} - \frac{\lambda \beta}{\Lambda(s + \lambda)} \right) N(s) - \frac{\Lambda}{a}\left(n(0) + \frac{\lambda c(0)}{s + \lambda}\right) = 0 $$

This is a linear, non-homogeneous, first order ODE of the form:

$$\frac{\mathrm{d}f(x)}{\mathrm{d}x} + P(x) f(x) - Q(x) = 0$$

Which has solutions as:

$$f(x) = \exp\left(-\int_0^{x} P(x') \mathrm{d}x\right) \left( \int_0^{x} \exp\left(\int_0^{x'} P(x'') \mathrm{d}x''\right) Q(x') \mathrm{d}x' + C \right)$$

So we can plug in our values:

$$N(s) = \exp\left(-\int_0^{s} \frac{\Lambda}{a}\left(s' + \frac{\beta}{\Lambda} - \frac{\lambda \beta}{\Lambda(s' + \lambda)} \right) \mathrm{d}s'\right) \left( \int_0^{s} \exp\left(\int_0^{s'} \frac{\Lambda}{a}\left(s'' + \frac{\beta}{\Lambda} - \frac{\lambda \beta}{\Lambda(s'' + \lambda)} \right) \mathrm{d}s''\right) \left(\frac{\Lambda}{a}\left(n(0) + \frac{\lambda c(0)}{s' + \lambda}\right)\right) \mathrm{d}s' + C \right)$$

But here we have a bit of a problem because we wind up with terms of $e^{s^2}$ under an integral. However, we still have an Inverse Laplace Transform to perform.

So my question is can this integral be solved? And if not, can we sidestep it with a clever application of the Inverse Laplace Transform?

Is there another way to approach this problem entirely that I'm overlooking?

For clarity, I'm seeking $n(t)$.

Answer 1

You can avoid having to use Laplace transforms altogether, and do the following:

1. Differentiate the $\frac{\text{d} n}{\text{d} t}$-equation to $t$,
2. Use the $\frac{\text{d} c}{\text{d} t}$-equation to express $\frac{\text{d} c}{\text{d} t}$ in terms of $c(t)$ and $n(t)$ in the result of step 1,
3. Use the $\frac{\text{d} n}{\text{d} t}$-equation to express $c(t)$ in terms of $n(t)$ and $\frac{\text{d} n}{\text{d} t}$ in the result of step 2,
4. Solve the resulting ODE for $n(t)$ in terms of confluent hypergeometric functions; see here or here for more information on these special functions.

To see that the resulting ODE is indeed of the required 'confluent hypergeometric form', first write $n(t) = e^{\lambda t} \,\nu(t)$. Then, in the resulting ODE for $\nu$, use the new variable $$ \tau = \frac{a t - \beta - \lambda \Lambda}{\sqrt{2 a \Lambda}}. $$ The result is an expression for $n(t)$ of the form $$ n(t) = e^{\lambda t} \left[c_0\,_1F_1\left(\frac{1}{2}-\frac{\beta \lambda}{2a},\frac{1}{2},\tau^2\right) + \tau\,c_1\,_1F_1 \left(1-\frac{\beta \lambda}{2a},\frac{3}{2},\tau^2\right) \right]. $$

Answer 2

I would try this: first you have

$$c(t)=e^{\lambda t} \left ( c(0) + \int_0^t e^{-\lambda s} \frac{\beta}{\Lambda} n(s) ds \right ).$$

Substituting that into the equation for $n$ gives an integrodifferential equation:

$$\frac{dn}{dt}=\frac{at-\beta}{\Lambda} n(t)-\lambda e^{\lambda t} c(0) - \lambda \int_0^t e^{\lambda (t-s)} \frac{\beta}{\Lambda} n(s) ds.$$

Now you can take the Laplace transform, and the Laplace transform of the convolution integral will just be a product. So you get

$$sN-n(0)=-\frac{a}{\Lambda} \frac{dN}{ds} - \frac{\beta}{\Lambda} N - \lambda \frac{c(0)}{s-\lambda} - \frac{\lambda \beta}{\Lambda} \frac{1}{s-\lambda} N$$ which in the usual form is $$ \frac{dN}{ds} + \frac{\Lambda}{a} \left ( s +\frac{\beta}{\Lambda} + \frac{\lambda \beta}{\Lambda} \frac{1}{s-\lambda} \right ) N = \frac{\Lambda}{a} \left ( -n(0) + \lambda \frac{c(0)}{s-\lambda} \right ).$$

At least in principle that's not hard to solve for $N(s)$ but I don't know how easy it will be to invert.

Incidentally, I would strongly encourage you to stick to definite integrals as I did here.