# 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 :)

## Answers and Replies

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.