Software for crank-nicholson method

  • Thread starter NiclasN
  • Start date
  • #1
4
0

Main Question or Discussion Point

Hi,

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

[tex]-\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)[/tex]

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

  • #2
4
0
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"?
 
  • #3
649
2
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.
 
  • #4
4
0
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.
 
  • #5
649
2
I think that in the Matlab PDE toolbox you must specify your equation in the form

[tex] - \nabla ( c \nabla N) + a N = f [/tex]

(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!
 
  • #6
4
0
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.
 

Related Threads on Software for crank-nicholson method

Replies
8
Views
2K
  • Last Post
Replies
0
Views
2K
Replies
4
Views
4K
Replies
23
Views
795
  • Last Post
Replies
14
Views
2K
  • Last Post
Replies
3
Views
4K
  • Last Post
Replies
4
Views
6K
Replies
0
Views
6K
Replies
1
Views
4K
Top