How to Solve a Highly Non-Linear PDE in Matlab

Hi all,

coming from the theoretical physics community I am not very much familiar with the numerical PDE and ODE solving capabilities of MATLAB.

I have a heavily non-linear elliptic PDE of the following form:

In the left hand side, it's the Laplacian in spherical coordinates; I am also interested in solving this equation in the case of beyond spherical symmetry.

In the right hand side, the prefactor outside the parenthesis is a constant number, n>=1, rho_c is another constant, rho(x) is a given function of x.

I need to impose boundary conditions in the center (the first derivative = 0), and very far away from the center (phi = some given number).

I have three questions:

1. how realistic is it for MATLAB to solve the given equation in 1D? Which function is it the best to use.

2. If I have a stiff system, namely rho_c is very-very much larger than the typical value of rho, is there any chance that MATLAB can handle this?

3. Moving to more complicated geometries, is this a difficult problem for methods like FEM?

I mainly want to know if the given equation is an easy one for engineering standards and whether it can be solved with MATLAB.

Thanks in advance.

Answer 1

[Precise Simulation]

You should fairly easily be able to enter this into the FEATool Multiphysics FEM toolbox as a custom PDE , for example the following code

fea.sdim  = { 'x' };
fea.grid  = linegrid( 10, 0, 1 );
fea.dvar  = { 'phi' };
fea.sfun  = { 'sflag1' };
fea.eqn   = parseeqn( 'phix_x + 2/x*phix_t = mc^2*Rc^2/(n+1)*(rho/rhoc-1/phi^(n+1))', fea.dvar, fea.sdim );
fea.bdr.d = { []  1 };
fea.bdr.n = {  0 [] };
fea.const = { 'rho'  '1+x' ;
        'rhoc' 1 ;
        'n'    1 ;
        'mc'   1 ;
        'Rc'   1 };
fea = parseprob( fea );
fea.sol.u = solvestat( fea, 'init', {1} );
postplot( fea, 'surfexpr', 'phi' )

should set up your problem with arbitrary test coefficients. Whether your actual problem is too nonlinear to converge is another issue though.

Answer 2

[Torsten]

Try "pdepe" with m=2.

It uses the method-of-lines approach based on a FEM discretization with a stiff ODE integrator in the background.

Best wishes

Torsten.