I'm trying to solve this problem numerically in Matlab:
$ \left\{ \begin{array}{rl} \frac{\partial P}{\partial t} &= \frac{\partial^2 P}{\partial x^2} \ \ \ (\star) \\ P(x,0) &= 1 \\ \frac{\partial P(0,t)}{\partial x} &= kP \\ P(1,t) &= 1 \end{array} \right. $
I considered a discretization of the space $[0,1]$ in $M$ intervals and time $[0,1]$ in $N$ intervals.
First I built a matrix $P$ to save the initial conditions, then for $P(x,0)=1$, I made: $P(m,0) = 1$, and for $P(1,t)$, I made $P(1,n) = 1$.
I used the $\theta$-method:
\begin{align*} \frac{\partial P}{\partial t} &\approx \frac{P^{n+1}_{m}-P^{n}_{m}}{\Delta t}\\ \frac{\partial^2 P}{\partial x^2} & \approx \frac{\theta P^{n+1}_{m+1}+(1-\theta)P^{n}_{m+1} -2(\theta P^{n+1}_{m}+(1-\theta)P^{n}_{m})+(\theta P^{n+1}_{m-1}+(1-\theta)P^{n}_{m-1}) }{(\Delta x)^2} \end{align*}
to obtain the Scheme:
\begin{equation} \begin{array}{l} \frac{-\theta}{(\Delta x)^2}P^{n+1}_{m-1}+\left( \frac{1}{\Delta t}+\frac{2\theta}{(\Delta x)^2} \right)P^{n+1}_{m}- \frac{-\theta}{(\Delta x)^2} P^{n+1}_{m+1} \\ =\frac{1-\theta}{(\Delta x)^2}P^{n}_{m-1}+\left( \frac{1}{\Delta t}-\frac{2(1-\theta)}{(\Delta x)^2} \right)P^{n}_{m}+ \frac{1-\theta}{(\Delta x)^2}P^{n}_{m+1} \end{array} \end{equation}
This scheme builds a system $Ax=b$, where $A$ (tridiagonal) has size $(M+1)\times (M+1)$ and for the first step with fixed time $n=0$, is posible to obtain the unknowns $P^1_{m-1}$, $P^1_{m}$, $P^1_{m+1}$, because by the inicial conditions, $P^0_{m-1}$, $P^0_{m}$ and $P^0_{m+1}$ are given.
But, I used approximations of order 2, so for the Neumann condition, I need to use an aproximation of order 2: "mirror imaging technique", i.e: If $J(P)$ is the flux: \begin{align*} J(P)_0 &= kP_0 \\ J(P)_0 &= \frac{J(P)_{-\frac{1}{2}}+J(P)_{\frac{1}{2}}}{2} \\ J(P)_{i+\frac{1}{2}} &= \frac{P_{i+1}-P_i}{\Delta x} \end{align*} with $i=0$: \begin{align*} J(P)_0 &= kP_0 \\ J(P)_{-\frac{1}{2}} &= 2J(P)_0 + J(P)_{\frac{1}{2}} \\ J(P)_{\frac{1}{2}} &= \frac{P_{1}-P_0}{\Delta x} \end{align*} Substituting in equation $(\star)$: \begin{align*} \frac{P^{n+1}_0-P^n_0}{\Delta t} -\frac{J(P)_{\frac{1}{2}}-J(P)_{-\frac{1}{2}}}{\Delta x} &= 0\\ \frac{P^{n+1}_0-P^n_0}{\Delta t} -\frac{J(P)_{\frac{1}{2}}-(2J(P)_0-J(P)_{\frac{1}{2}})}{\Delta x} &= 0\\ \frac{P^{n+1}_0-P^n_0}{\Delta t} -\frac{\frac{P_1-P_0}{h}-(2kP_0-\frac{P_1-P_0}{\Delta x})}{\Delta x} &= 0\\ \frac{P^{n+1}_0-P^n_0}{\Delta t} -\frac{2kP_0}{h} - 2\frac{(P_1-P_0)}{\Delta x^2} &= 0\\ \frac{P^{n+1}_0-P^n_0}{\Delta t} -\frac{2kP_0}{h} - 2\frac{P_1}{\Delta x^2} + 2\frac{P_0}{\Delta x^2} &= 0 \\ \frac{P^{n+1}_0-P^n_0}{\Delta t} -2\frac{P_1}{\Delta x^2} - \frac{2}{\Delta x}(k+\frac{1}{h})P_0 &= 0 \\ \end{align*} Replacing expressions of $\theta$-method in the above equation: $$ \frac{P^{n+1}_0-P^n_0}{\Delta t} -\frac{2}{\Delta x^2}\left( \theta P^{n+1}_1 + (1-\theta)P^n_1 \right) - \frac{2(k+\frac{1}{h})}{\Delta x}\left( \theta P^{n+1}_0 + (1-\theta)P^n_0 \right) = 0 $$ Then, grouping terms: $$\left( \frac{1}{\Delta t} + \frac{2(k+\frac{1}{\Delta x})\theta}{\Delta x}\right) P^{n+1}_0 - \frac{2\theta }{h^2} P^{n+1}_1 = \left( \frac{1}{\Delta t} - \frac{2(k+\frac{1}{h})(1-\theta)}{\Delta x} \right) P^n_0 + \frac{2(1-\theta)}{\Delta x^2} P^n_1 $$
Finally using the above equation, the scheme and the another initial condition, I fill entries of $P$ (except those corresponding to the boundary conditions) with row vectors that are solutions of the problem for a fixed time.
Here my questions:
- Is the deduction for the Neumann condition well done?
- How does it affect my program in matlab the Neumann condition? I have particular problems with that because I first thought that only affected the first equation of the system $ Ax = b $, but in that case, the system would have more unknowns than equations and could not solve the system.
Thanks for your help.