Solve heat equation with nonlinear term in matlab - Does my solution look reasonable?

Question

4) Write a program to solve the following initial boundary value problem: $$ v_t = 0.01 v_{xx} + 100 v (v-1) (.25 - v) $$ with $v_x(0,t) = v_x(1,t) = 0$ for all $t>0$, and $$ v(x,0) = \left\lbrace\begin{aligned} 0\quad &\text{if } 0 \leq x \leq .40 \\ a\quad &\text{if } .40 < x < .60 \\ 0\quad &\text{if } .60 \leq x \leq 1.0 \end{aligned}\right. $$ Solve this problem with $a = .20$ and with $a = .90$ for the time interval $0\leq t \leq 1$ and comment on the difference in the nature of the solution. Try to explain this difference in behavior.

I know how to solve $v_t=Cv_{xx}$ and I think I am ok on the Neumann BC ad the piece wise component, but I don't know how to account for the nonlinear term though. I imagine that I might want to use newton's method or Runge Kutta methods or something but I don't know what I would want to solve for. I am really just at a loss on this problem :(

UPDATE: I have found that to satisfy the Neumann conditions we need the following matrix I think:

Any idea if this looks reasonable for $a=0.3$? Each line is at an updated time.

Answer 1

Let us use an explicit FD method with centered differentiation in space, where $k$ is the time-step and $h$ the mesh size. Thus, if $v_j^n \simeq v(j h, n k)$, we have $$ v_j^{n+1} = (1 - 2\lambda) v_j^{n} + \lambda (v_{j+1}^{n} + v_{j-1}^{n}) + k R(v_j^{n})\, , $$ where $R(v) = 100 v (v-1) (.25 - v)$ is displayed below, and $\lambda = 0.01 k/h^2$ is the Fourier number. The boundary conditions are discretized using centered finite differences: \begin{aligned} v_{-1}^n &= v_{+1}^n \\ v_{J_x+1}^n &= v_{J_x-1}^n \end{aligned} where $J_x = 1.0/h$.

In the case $R(v) = 0$, the numerical method is stable and convergent under the condition $\lambda \leq 1/2$. In the nonlinear case, the method writes as a nonlinear dynamical system ${\bf v}^{n+1} = {\bf f}({\bf v}^{n})$ with ${\bf f}$ nonlinear, which makes stability analysis more difficult.

Here is a Matlab code which implements the method, and the numerical solution for $a = .20$ and $a = .90$ at the final time $t = 1.0$:

% relaxation function
R = @(u) 100*u*(u-1)*(0.25-u);

% domain
Jx = 200;

% initial data
a = 0.9;

% data
h = 1/Jx;
x = -h:h:1+h;
lambda = 0.25;
k = h^2*lambda/0.01;
v = zeros(size(x));
v((0.4<x)&(x<0.6)) = a;
vtemp = v;

% initialization
figure;
clf;
h = plot(x,v);
xlim([0 1]);
ylim([0 1]);
xlabel('x');
ylabel('v');
t = 0;
ht = title(sprintf('t = %0.3f',t));

% loop
while t+k<1
    % scheme
    for j=2:Jx+2
        vtemp(j) = (1-2*lambda)*v(j) + lambda*(v(j-1)+v(j+1)) + k*R(v(j));
    end
    % boundary conditions
    vtemp(1)    = vtemp(3);
    vtemp(Jx+2) = vtemp(Jx);

    t = t+k;
    v = vtemp;
    set(h,'YData',v);
    set(ht,'String',sprintf('t = %0.3f',t));
    drawnow;
end

Let us examine the relaxation part $v_t = R(v)$ (cf. splitting). From the graph of $R$, one can deduce that the long-time solution to the IVP $v(0) = .20$ is $v(\infty) = 0$, while the long-time solution to the IVP $v(0) = .90$ is $v(\infty) = 1.0$. This explains the different nature of the solutions for $a= .20$ and $a = .90$.