Software for crank-nicholson method

[NiclasN]

Hi,

I try to numerically solve the following partial differential equation for N(r,z) with a Dirichlet boundary condition.

$$-\frac{\partial^2N}{\partial r^2}-\frac{\partial^2N}{\partial z^2} + f(r,z) \frac{\partial N}{\partial r} + g(r,z) \frac{\partial N}{\partial z} = h(r,z)$$

Mathematica is apparently not able to do it, because it is not an initial value problem. I also wasn't successful with the MATLAB pde tool. The Crank-Nicolson approximation seems to be the right way to go.

My question is which is the best software for solving this problem, so that I don't have to implement the algorithm myself.

Thanks a lot :)

 

[NiclasN]

Maybe I should have colled the thread more accurately "numerical solution to a 2nd order pde with dirchlet boundary conditions". Is there no software that can solve this without discretizing it "by hand"?

 

[torquil]

By "Matlab PDE tool", do you mean the Matlab PDE toolbox? How did it fail? I was under the impression that it should be able to solve such a problem.

 

[NiclasN]

Yes, I meant the toolbox. I failed because I could not specify this form of equation. It's not a time dependent equation in the sense that you have one variable (t) that occurs only in first derivative. There is first and second derivatives in both variables. Also I think for the toolbox you have to write the equation in "gradient form" for cartesian variables.
The problem stems from a steady state diffusion equation in cylindrical coordinates, symmetric in phi. The boundary condition is N=0 on the boundary of the cylinder.
I would of course be glad if I was wrong and it can be done with Matlab.

 

[torquil]

I think that in the Matlab PDE toolbox you must specify your equation in the form

$$- \nabla ( c \nabla N) + a N = f$$

(where f is not the same as your f). I think you can find c such that your equation can be written on that form. And the matlab-f would be c*h where h is your h... The first order derivatives on your N would come from the term where the nabla above acts on c, which would determine c.

The (r,z) domain is a rectangle. You know the Dirichlet boundary condition N=0 e.g. at r=1 and z=0 and z=1, since that is the boundary of the solid cylinder you are working with (modulo units of length...)

You also need a condition at r=0. But by your angular symmetry, isn't that just a Neumann condition at r=0, since the r-derivative of N is zero at r=0 since the 3d problem is rotationally invariant about the z-axis?

Its just something to consider, if you hadn't before, or if you hadn't seen that form of the equation to be specified in Matlab. Good luck!

 

[NiclasN]

Thanks a lot! I did not consider that c can be a matrix with coefficients dependent on (r,z) i.e. (x,y) but apparently this is possible so maybe you are right. I'll try again.