Solving Partial Differential Equation using Crank-Nicolson Method in MATLAB

[electronic engineer]

Hello,

i'm willling to solve partial differential equation of the form:

dy/dt = d2y/dx2 + 2/x * dy/dx - 4y

where d here means partial derivation of course.

boundary conditions are: dy/dx(0,t) = 0, y(1,t) = 1, t>0

initial condition is: y(x,0) = 0;

(x,y) is position and t is time.

and t should be also chosen so as the solution doesn't differ from the stationary solution more than 5%, where the stationary solution is the one who is independent on time t.

Can you please give me some hints on how to start or think about this problem? I need to build some code to solve this equation using matlab.

Thank you. I'm looking forward to your answers.

 

[bigfooted]

Hi,

Have you already programmed the Crank-Nicolson method in matlab? You can then play around with it and get a feeling for what's going on and how the stepsize changes the long-term solution. You could post the code here if you have problems getting it running, it should be like 20 lines or so, but please also add comment lines if you post it.
Have you already found the exact solution of the steady-state problem? You need this to compare your numerical solution at t=large.

 

[electronic engineer]

Hi, thanks for your reply. I haven't programmed it in MATLAB yet. But I'm not familiar how to solve it generally and how to hav the final formula for the solution so then i would start programming. The main is how to get it resolved. That's the point.

 

[bigfooted]

You could start with the analytic solution of the steady-state equation. I guess you followed some 'Intro to ODE's' course, so you can first put dy/dt=0 and solve the linear second order ODE. The theory is explained in books like the one from Boyce and DiPrima, 'Elementary Differential Equations...'

Do you have any experience with MATLAB or numerical methods? If not, I would start with a simpler problem,
y'(x)=y(x), y(0)=1 on x=[0,1]

Solve the problem analytically and then solve it numerically using the forward Euler method: y(x+1)=y(x)+dx*y'(x)
http://en.wikipedia.org/wiki/Euler_method , see the example

y(0)=1, dy(0)=y(0);
y(0.01) = y(0)+0.01*dy(0)

and check the numerical solution with the exact solution.

The Crank-Nicolson method is a combination of the forward and backward Euler method:
http://en.wikipedia.org/wiki/Crank–Nicolson_method

 

[electronic engineer]

Thank you for the reply Bigfooted.

I'm familiar with the numerical methods. I have applied the crank-nicolson relation so after many corrections I had some formula. But still I'm confused about how to move forward. Or maybe I don't remember verywell when it comes to more difficult tasks. I'll keep trying!

 

[Chestermiller]

Thank you for the reply Bigfooted.

I'm familiar with the numerical methods. I have applied the crank-nicolson relation so after many corrections I had some formula. But still I'm confused about how to move forward. Or maybe I don't remember verywell when it comes to more difficult tasks. I'll keep trying!

You are aware that you need to discretize the differential equation into finite differences in the spatial direction to obtain a set of coupled first order ODEs, correct?

 

[electronic engineer]

You are aware that you need to discretize the differential equation into finite differences in the spatial direction to obtain a set of coupled first order ODEs, correct?

yes, so what are you aiming to?

 

[Chestermiller]

yes, so what are you aiming to?

If this is the case, it isn't clear to me what your difficulty is in applying Crank-Nicholson to solve a set of coupled first order linear ordinary differential equations. Can you articulate what you are having difficulty with?

Chet

 

[electronic engineer]

If this is the case, it isn't clear to me what your difficulty is in applying Crank-Nicholson to solve a set of coupled first order linear ordinary differential equations. Can you articulate what you are having difficulty with?

Chet

I will tell you, I applied Crank-Nicholson directly to solve the second order..as some examples i read in some textbook but the examples were solved there in a completely different way that's way i had problem.. I'm familiar with the idea that you said but not sure how to apply it to this concrete example!

 

[Chestermiller]

I will tell you, I applied Crank-Nicholson directly to solve the second order..as some examples i read in some textbook but the examples were solved there in a completely different way that's way i had problem.. I'm familiar with the idea that you said but not sure how to apply it to this concrete example!

Was the example in your book a single ode or was it also a set of coupled odes?

Chet

 

[electronic engineer]

Was the example in your book a single ode or was it also a set of coupled odes?

Chet

the example was a single ode!

 

[Chestermiller]

the example was a single ode!

So then, is your question, "How do I apply this to a set of coupled odes?"

It's very hard to help you if you don't reveal your difficulty in sufficient detail.

Chet

 

[electronic engineer]

So then, is your question, "How do I apply this to a set of coupled odes?"

It's very hard to help you if you don't reveal your difficulty in sufficient detail.

Chet

I would prefer to use the set of coupled odes as was in the example of my textbook. But it can also be any method that make the solution easier so i don't really know! If it's less difficult to solve it by single ode then i would go for it also. I'm not sure which one can go work better for this example.

You are of course right. I'm trying to study many materials in order to be more oriented in this. I was following some way but then a block stood in the way.

 

[electronic engineer]

I think this script can work very well for my example:

http://www.physycom.unibo.it/bazzani_web_dir/modelli_numerici/algoritmipde.pdf

 

[Chestermiller]

I think this script can work very well for my example:

http://www.physycom.unibo.it/bazzani_web_dir/modelli_numerici/algoritmipde.pdf

Yes. This write up looks very good. Be very careful how you handle the boundary condition at x = 0.

Chet