PDE Matlab Help

Spitfir3 · 17 Jun 2015 at 19:24

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

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

skaif · 17 Jun 2015 at 19:39

Better trying this in the programming section. GD will only get you pointless answer's

!

Scotty123 · 17 Jun 2015 at 19:48

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?

Spitfir3 · 18 Jun 2015 at 11:31

Scotty123 said:
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?

Scotty123 · 18 Jun 2015 at 14:02

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.

Spitfir3 · 18 Jun 2015 at 16:04

Scotty123 said:
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.

Scotty123 · 18 Jun 2015 at 16:41

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).