Crank-Nicolson method for a parabolic differential equation

[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

 

[Astronuc]

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

See if these help

http://math.fullerton.edu/mathews/n2003/CrankNicolsonMod.html

Nice summary of PDE's
instruct1.cit.cornell.edu/courses/ engr241/LectNotes/9_PDEs_Fa00_notes.doc use save target as

 

[jamesponty]

Thanks for your reply. I've looked through the links you gave me fairly thoroughly but still really don't 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".

 

[LeBrad]

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

 

[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 can't see how to alter the time interval. It is also not reaching a steady state either. I have just added the -2 as constant because it gives a form in the textbook 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

 

[jamesponty]

Can anyone offer a suggestion please?

 

[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.

 

[LeBrad]

Crank Nicholson says something like

<br /> u^{j+1} = \frac{1}{2}[f(u^j(x),t^{j})+f(u^{j+1}(x),t^{j+1})]<br />

And you are solving

<br /> \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} - 2<br />

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

<br /> \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<br />

Then rearrange

<br /> 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} <br /> + <br /> \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}<br /> - 2 \Delta t<br />


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.

 

[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