# Crank-Nicolson method for a parabolic differential equation

1. Sep 10, 2005

### jamesponty

Hi,

Just wondering if somebody could point me in the right direction.

I have the parabolic differential equation:

du/dt=d2u/dx2 -2

I need to get a MATLAB script to solve this but I have no idea how to handle the -2. I also have no idea how to edit the MATLAB scripts I've found to encorporate this???

Any ideas?

Thanks
Jamesponty

2. Sep 10, 2005

### Staff: Mentor

Last edited by a moderator: May 2, 2017
3. Sep 13, 2005

### jamesponty

Thanks for your reply. I've looked through the links you gave me fairly thoroughly but still really dont know where to start from. I have sample code for Matlab to solve the problem du/dt=d^2u/dx^2 and am only meant to be able to alter this code to account for the "-2".

4. Sep 13, 2005

Just add -2 at ever iteration. Why don't you post what you have so far.

5. Sep 17, 2005

### jamesponty

OK, here is what I've done so far. It is only an adaption of code that we were given for the assignment. We are supposed to plot a number of representitive prolfiles for temp and use trial and error to determine a suitable value for the t interval. It also states to use N=12 steps and to advance the solution until you consider steady state to be reached.

I am pretty sure I have adjusted the code to use 12 steps but cant see how to alter the time interval. It is also not reaching a steady state either. I have just added the -2 as constant becuase it gives a form in the text book of A*u(j+1)=Bu(j)+c. I have just made c = -2.

The other relevant info is the boundary conditions
u(x,0)=100x(1-x)^2 for 0<x<1
u(0,t)=u(1,t)=0 for t>0

Where should I go from here? Am I on the right track?

hold on
r=0.7;
cons=-2;
A=eye(11);
A=2*(1+r)*A;
for i=1:10
A(i,i+1)=-r;
A(i+1,i)=-r;
end
B=eye(11);
B=2*(1-r)*B;
for i=1:10
B(i,i+1)=r;
B(i+1,i)=r;
end
C=inv(A)*B+inv(A)*cons;
y=0:1/12:1.001;
x=y(2:12);
u=(100*x.*(1-x).^2)';
v(1)=0;
v(13)=0;
for j=1:31
for i=1:11
v(i+1)=u(i);
end
plot(y,v)
u=C*u;
end

Last edited: Sep 18, 2005
6. Sep 22, 2005

### jamesponty

Can anyone offer a suggestion please?

7. Sep 22, 2005

### saltydog

What are the initial conditions and boundary conditions? Can you set it up like this please and fill in any changes that are required?

$$\text{DE:}\quad u_{t}=u_{xx}-2 \quad 0\leq x \leq L$$

$$\text{BC:}\quad u(0,t)=0\quad u(L,t)=0$$

$$\text{IC:}\quad u(x,0)=f(x)$$

Me, first I'd solve it directly using Duhamel's method to get the right answer, then debug the Crank-Nicoloson method to match what I got via separation of variables. Hey, I ain't proud.

Edit: Oh yea, if for some reason your problem is such that it's not solved easily directly via SOV and Duhamel's method, then set up one that IS easily solved directly, get the right answer, then debug Crank-Nicoloson to match that one, then substitute your problem into the bug-free method.

Last edited: Sep 22, 2005
8. Sep 22, 2005

Crank Nicholson says something like

$$u^{j+1} = \frac{1}{2}[f(u^j(x),t^{j})+f(u^{j+1}(x),t^{j+1})]$$

And you are solving

$$\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} - 2$$

So discretize your equation, I use subscript i for position and superscript j for time.

$$\frac{u_i^{j+1}-u_i^j}{\Delta t} = \frac{u_{i+1}^{j}-2u_i^j+u_{i-1}^j}{2 \Delta x^2} + \frac{u_{i+1}^{j+1}-2u_i^{j+1}+u_{i-1}^{j+1}}{2 \Delta x^2} -2$$

Then rearrange

$$u^{j+1}_i = \frac{\Delta t}{2 \Delta x^2}u^{j}_{i+1} + (-2 \frac{\Delta t}{2 \Delta x^2}+1)u^{j}_{i} + \frac{\Delta t}{2 \Delta x^2}u^{j}_{i-1} + \frac{\Delta t}{2 \Delta x^2}u^{j+1}_{i+1} + (-2 \frac{\Delta t}{2 \Delta x^2})u^{j+1}_{i} + \frac{\Delta t}{2 \Delta x^2}u^{j+1}_{i-1} - 2 \Delta t$$

That is an equation for an interior node in your domain.

I'm pretty sure I see [one of] the problem with your code. I guess just telling you to add -2 at each step is not entirely accurate. It is weighted by something, and you don't want it on every node because of boundary conditions. See if you can figure out what to do from there.

9. Dec 17, 2009

### dumurum

hi,
i'm stuck with that homework .could someone help me pls?

d2u/dt2=d2u/dx2 , 0<x<1 , t>=0,

boundary conditions : u(0,t)=0=u(1,t), t>=0
initial conditions : u(x,0)=sin(pi*x) , 0<x<1
u(x,0)=0, 0<x<1