Using MATLAB to solve a system of two PDEs encountered in Mathematical Biology

Question

I am in an undergraduate summer project and am studying a paper, I hope I could solve the following two PDEs with two boundary conditions for $y_1(t,a)$ and $y_2(t,a)$ using MATLAB. $$\frac{\partial y_1}{\partial t}+\frac{\partial y_1}{\partial a}=f_1(t)y_1, \qquad t\geqslant0, \quad a\in[0,\tau_1],$$ with boundary condition $y_1(t,0)=c_1$. $$\frac{\partial y_2}{\partial t}+\frac{\partial y_2}{\partial a}=f_2(t)y_2, \qquad t\geqslant0,\quad a\in[0,\tau_2],$$ with boundary condition $y_2(t,0)=c_2y_1(t,\tau_1)$.

Here, $y_1(t,a)$ and $y_2(t,a)$ are functions of $t$ and $a$. The functions $f_1(t)$ and $f_2(t)$ are known functions of $t$ but are rather complicated, and $\tau_1,\tau_2,c_1,c_2$ are known constants.

Since $y_1(t,a)$ is involved in the boundary condition of the second PDE for $y_2(t,a)$, I aim to solve the first PDE for $y_1(t,a)$ first, and then use the results in solving the second PDE for $y_2(t,a)$.

I am new in both PDEs and MATLAB, although I have spent several days searching for possible ways in doing it, I still have got no hints. It would be greatly appreciated if some hints or possible directions I could look into could be provided. Thanks in advance.

Answer 1

First order PDEs can generally be solved by the method of characteristics. To solve your first PDE using this method, put $$ \phi_z(x) = y_1(z + x, x)\ . $$ Then \begin{eqnarray} \frac{d\phi_z}{dx}(x)&=&\frac{\partial y_1}{\partial t}(z+x,x)+\frac{\partial y_1}{\partial a}(z+x,x)\\ &=& f_1(z+x)\,\phi_z(x)\ . \end{eqnarray} This is a first-order ODE for $\ \phi_z(x)\ $, which has the solution: \begin{eqnarray} \frac{\phi_z(x)}{\phi_z(0)}&=&e^{\int_0^xf_1(z +u)\,du}\\ &=& e^{\int_z^{x+z}f_1(u)\,du}\ . \end{eqnarray}

Now, $\ \phi_z(0)=y_1(z, 0) = c_1\ $, and $\ y_1(t,a)=\phi_{t-a}(a)\ $, so $$ y_1(t,a)=c_1 e^{\int_{t-a}^tf_1(u)\,du}\ . $$ Your second PDE can now be solved similarly. When solving the corresponding ODE, however, you will need to integrate over the interval from $\ x\ $ to $\ \tau_1\ $, rather than from $\ 0\ $ to $\ x\ $ as I did above.

Answer 2

The method of characteristics provides the solution in the domain $t\geq a \geq 0$. The characteristic curves along which the information $\frac{\text d}{\text d t} y_k = f_k y_k$ propagates are the straight parallel lines $x = t-t_0$ with $t_0 \geq 0$. Hence, we have $$ y_1(t,a) = c_1 \exp\left(\int_{t-a}^t f_1(s)\, \text d s\right) = c_1 \exp\left(\int_{0}^a f_1(t-\tau)\, \text d \tau\right) . $$ Similarly, we have for $t\geq a \geq 0$ $$ y_2(t,a) = c_2\, y_1(t,\tau_1) \exp\left(\int_{0}^a f_2(t-\tau)\, \text d \tau\right) $$ where $y_1(t,\tau_1)$ with $t\geq \tau_1$ can be deduced from the previous step. Of course, this solution can then be implemented in Matlab, e.g. by using the Matlab function integral for the numerical integration of $f_k$. Alternatively, an iterative numerical resolution by finite-volume methods can be implemented.