Matlab: Solving $u_t = f(u) u_{xx} + g(u) u_x + h(u) u$

Question

I am trying to solve the following equation numerically:

$$\dfrac{\partial u}{\partial t}=\dfrac 3 x\dfrac \partial{\partial x}\left(\sqrt x \frac \partial {\partial x}(\nu(u,u^2,\cdots,u^k)u\sqrt x)\right)$$

Expressed in its expanded form, this is:

$$u_t=3\nu(u,u^2,\cdots)u_{xx}+\left(6\nu_x(u,u^2,\cdots)+\dfrac{3\nu(u,u^2,\cdots)}{2x}\right)u_x+\left(\nu_x(u,u^2,\cdots)\left(3+\dfrac 1 {2x}\right)+\dfrac{3\nu(u,u^2,\cdots)}{x}\right)$$

Where $u$ is a function of $u=u(t,x)$. $\ \nu$ a is a nonlinear function of $u$ and another variable, but for a given $u$ I can obtain it.

What is the best way in terms of stability, performance, and accuracy to solve this? I was trying method of lines in MATLAB, but I have not been able to make any progress.

Answer 1

If you know for sure that your PDE is well posed, then I would suggest the following:

1. Try discretizing using an explicit finite difference scheme first.

2. If you find that you get unreasonable / unstable solutions, try to make your temporal stepsizes very small in comparison to your spatial stepsizes (by several orders of magnitude).

3. If this still doesn't work, try using an implicit finite difference scheme. Since your equations are non-linear, it would be best to solve using newton's method.

4. If all else fails, check that your problem is well-posed. You may want to try proving/disproving that solutions to your problem exist and are unique for the initial/boundary conditions you impose.

Answer 2

It's hard to give concrete advice without knowing how your $\nu(u,u^2,\dots)$ really looks, but if you're in Matlab, you may want to try the Chebfun system, which uses collocation spectral methods and automatic differentiation to solve linear and non-linear ODEs and, in some cases, PDEs.

If you install Chebfun, the fastest way to try your example is to start the chebgui, a graphical user interface for ODEs and PDEs, and look at some of the PDE examples and modify them for your own function.

Disclaimer: I am part of the Chebfun developer team :)

Update

Following Aron Ahmadia's request, I've cooked-up a small example showing how to compute $\partial u/\partial t$ for $\nu(u) = u + u^2$ with boundary conditions $u(1)=u(2)=1$:

% Make some (any) nu(u)
nu = @(u) u + u.^2;

% Set-up the differential equation for dudt
dudt = @(x,u) 3./x.*diff( sqrt(x) .* diff( nu(u).*u.*sqrt(x) ) );

% Create a chebop with that operator on [1,2]
N = chebop( dudt , [1,2] );

% Set some arbitrary boundary conditions
N.lbc = 1;
N.rbc = 1;

% Solve for some righ-hand side
u = N \ 1;
plot(u);

The syntax may be a bit cryptic, but, again, I'd recommend having a look at chebgui, which has a more intuitive syntax and is easier (and more fun) to play with.