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

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.

• What are your boundary conditions? – Paul Apr 2 '12 at 19:09
• What is the domain? – David Ketcheson Apr 3 '12 at 4:49
• Paul I want to solve it generally, but If you need some boundary condition can be this one: d/dx (vusqrt(x))=0 ) at x0 – Nikko Apr 3 '12 at 12:30
• David I am not sure what you mean. The spatial domain? from 0 to some finite arbitrary value. – Nikko Apr 3 '12 at 12:32

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.

• In step 1 above, keep the coefficients fixed, i.e. without the nonlinearity. In that case, you should be able to compute an exact solution and verify that your scheme works, before going to the nonlinear scheme. – Wolfgang Bangerth Apr 2 '12 at 19:39
• Ok. Well I have no problem with the coef are fixed actually. The problem is when they are not fixed. In fact "v" is a really messy function of "u". – Nikko Apr 8 '12 at 17:35
• By the way, which kind of implicit finite difference scheme do you think is better? – Nikko Apr 10 '12 at 17:55

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.

• Actually, I'm interested in seeing the code for solving this in Chebfun. If @Nikko does not update his question, do you mind "coming up" with some initial conditions and a $\nu$ then showing what a solution might look like? – Aron Ahmadia Apr 3 '12 at 8:15
• Hey! thanks for the help. Sorry for the delay but I am working on this since one months and I have a lot of work hehe. So, I will check this Chebfun and I see what happened. But for the record, I don't have problems with "easy" prescription for "v" such us v=a+a^2 or something. My problem comes because of this. My "v" function depends on "u" and another variables "T" which depends on "u" also. so i have this v= f(u,T) but T=G(u). The problem is that I can't have a simplier formula for T in terms of "u". T actually is a complicated expression. – Nikko Apr 8 '12 at 17:42
• @Nikko: It doesn't (shouldn't) matter how complicated $\nu(u)$ is as long as it is differentiable. – Pedro Apr 10 '12 at 8:00
• That is true. It shouldn't. – Nikko Apr 10 '12 at 12:24
• Pedro, how you write nu if nu where an equation which I need to solve? nu=f(T,x) and T=g(T,u) then I need to solve g(T,u)=0 – Nikko Apr 10 '12 at 18:04