PDE Matlab Help

Associate
Joined
9 Oct 2004
Posts
1,376
Location
Paris
Hi all,

Was wondering there's any programmers / people good with MATLAB on here who could help me.. i'm rather desperate. A crate of beer would be given to anyone who could help :D

I've got a 2D non-linear elliptical partial differential equation in the form

ux - uxx + uyy = 0

I need to solve using a finite difference scheme, with various dirichlet and neumann boundary conditions.

Tried using PDETOOL in Matlab but it uses FEM which is inconsistent with the other numerical schemes I'm using.

Cheers :)
 
Associate
Joined
16 Aug 2010
Posts
1,373
Location
UK
What's the boundary conditions and domain size?

It's simple enough to solve without resorting to any Matlab functions, can just do it yourself with the appropriate difference scheme. Get your coupled equations for the node values to give you your matrix, then invert.

Uniform mesh I assume?
 
Last edited:
Associate
OP
Joined
9 Oct 2004
Posts
1,376
Location
Paris
What's the boundary conditions and domain size?

It's simple enough to solve without resorting to any Matlab functions, can just do it yourself with the appropriate difference scheme. Get your coupled equations for the node values to give you your matrix, then invert.

Uniform mesh I assume?

I need to do it in Matlab as I need the values for other PDEs i'm using. Also would be handy to quickly adjust the mesh resolution - and yep uniform mesh.

https://www.dropbox.com/s/f7jv4gsfoys4azf/pde.pdf?dl=0

This is the problem (it's a 2D diffusion flame model), looking at a 10 x 1 domain size. The left hand side is the inlet; the top and bottom parts are the inlets for the oxidiser (Z=0) and the middle is the inlet for the fuel (Z=1). The top and bottom are walls (hence the y BC).

Think you could help?
 
Associate
Joined
16 Aug 2010
Posts
1,373
Location
UK
I didn't say you couldn't do it in Matlab? I said you don't need to use the Matlab functions for PDEs that blackbox the solver, since it's simple enough to solve via constructing your matrix and inverting, which of course you do in Matlab.

Using your equation, write out the difference scheme in terms of i and j for x and y respectively, e.g. the central difference approximation to dZ/dx is (Z^(i+1,j) - Z^(i-1,j)) / delta. You will then get an expression for Z in terms of the nodal positions. They are coupled of course, but can be put into matrix form (with the boundary conditions incorporated). Then you just invert to solve for Z at each node.
 
Last edited:
Associate
OP
Joined
9 Oct 2004
Posts
1,376
Location
Paris
I didn't say you couldn't do it in Matlab? I said you don't need to use the Matlab functions for PDEs that blackbox the solver, since it's simple enough to solve via constructing your matrix and inverting, which of course you do in Matlab.

Using your equation, write out the difference scheme in terms of i and j for x and y respectively, e.g. the central difference approximation to dZ/dx is (Z^(i+1,j) - Z^(i-1,j)) / delta. You will then get an expression for Z in terms of the nodal positions. They are coupled of course, but can be put into matrix form (with the boundary conditions incorporated). Then you just invert to solve for Z at each node.

Thanks.

I'm trying to modify this script accommodate the extra first order derivative however i'm not getting anywhere:

http://www.mathworks.com/matlabcentral/fileexchange/38091-2d-laplace-equation

Is the method you're describing the method on page 14 of this document: http://www.mit.edu/~lululi/school/16.920_Numerical_Methods_for_PDEs/notes/__all__.pdf ?

Cheers.
 
Associate
Joined
16 Aug 2010
Posts
1,373
Location
UK
Yes it is on page 14.

You know your expressions for the derivatives in terms of the finite differences, stick those into your equation, simplify if needed so you have collected like terms involving Z. This should make it obvious as to what the matrix will be, and you can construct it.

It doesn't matter about the first order derivative, as when you include its central difference along with the second order x derivative, the nodal values at i+1,j and i-1,j will be collected. They will have coefficients depending upon the Peclet number.

For example, the nodal value at Z^(i+1,j) will have a coefficient of (Pe . delta - 1) / (Pe . delta^2), where I have simply combined the coefficient 1 / delta - 1 / (Pe . delta^2).
 
Last edited:
Back
Top Bottom