I'm working on a transient 2D heat equation model and am having a few

problems with the boundary conditions for my 2D plate.

here's my code:

%% IMPLICIT CRANK NICOLSON METHOD FOR 2D HEAT EQUATION%%

clc;

clear all;

% define the constants for the problem

M = 25; % number of time steps

L = 1; % length and width of plate

k = 0.02; % time step

N = 5; % number of elements in plate (NxN)

h = L/(N-1); % grid spacing

x = [0:h:L]; % vector of x values

y = [0:h:L]; % vector of y values

%

[xx yy] = meshgrid(x,y);

%

% set constants a and c

a = (1 + (2*k)/(h^2));

c = (-k/(2*(h^2)));

% set the matrix for H

A = a*eye(N) + c*diag(ones(1, N-1), 1) + c*diag(ones(1, N-1),

-1);

B = c*eye(N);

C = zeros(N);

%

% Create the matrix H, composed of sub-matrices A,B and C

Ap = eye(N);

Bp = diag(ones(1, N-1), 1) + diag(ones(1, N-1), -1);

Cp = ones(N) - Ap - Bp;

%

% Define the H matrix using KRON to "replace" 1's in Ap, Bp, and Cp

with copies of A, B, and C %

H = kron(Ap, A) + kron(Bp, B) + kron(Cp, C);

%

% % set up initial conditions for plate

a = 4; % the initial cond value for the plate

u = a*ones(N,N);

%

% set the boundary conditions for the plate

u(1,:) = 0; % b.c. for left hand end of plate

u(:,1) = 0; % b.c. for right hand end of plate

u(:,N) = 0; % b.c. for top hand end of plate

u(N,:) = 0; % b.c. for bottom hand end of plate

%

u = u(:);

for m = 1:M

u = H\u

axis manual; axis([0 L 0 L 0 a]);

drawnow

uu = reshape(u,N,N)

surf(xx,yy,uu), %colorbar

end

I think I need to put the b.c.'s inside the main loop for them to be

'remembered' throughout the iteration process. At the moment it

forgets them as soon as the iteration begins.

Any ideas would be most appreciated,

Geoff

problems with the boundary conditions for my 2D plate.

here's my code:

%% IMPLICIT CRANK NICOLSON METHOD FOR 2D HEAT EQUATION%%

clc;

clear all;

% define the constants for the problem

M = 25; % number of time steps

L = 1; % length and width of plate

k = 0.02; % time step

N = 5; % number of elements in plate (NxN)

h = L/(N-1); % grid spacing

x = [0:h:L]; % vector of x values

y = [0:h:L]; % vector of y values

%

[xx yy] = meshgrid(x,y);

%

% set constants a and c

a = (1 + (2*k)/(h^2));

c = (-k/(2*(h^2)));

% set the matrix for H

A = a*eye(N) + c*diag(ones(1, N-1), 1) + c*diag(ones(1, N-1),

-1);

B = c*eye(N);

C = zeros(N);

%

% Create the matrix H, composed of sub-matrices A,B and C

Ap = eye(N);

Bp = diag(ones(1, N-1), 1) + diag(ones(1, N-1), -1);

Cp = ones(N) - Ap - Bp;

%

% Define the H matrix using KRON to "replace" 1's in Ap, Bp, and Cp

with copies of A, B, and C %

H = kron(Ap, A) + kron(Bp, B) + kron(Cp, C);

%

% % set up initial conditions for plate

a = 4; % the initial cond value for the plate

u = a*ones(N,N);

%

% set the boundary conditions for the plate

u(1,:) = 0; % b.c. for left hand end of plate

u(:,1) = 0; % b.c. for right hand end of plate

u(:,N) = 0; % b.c. for top hand end of plate

u(N,:) = 0; % b.c. for bottom hand end of plate

%

u = u(:);

for m = 1:M

u = H\u

axis manual; axis([0 L 0 L 0 a]);

drawnow

uu = reshape(u,N,N)

surf(xx,yy,uu), %colorbar

end

I think I need to put the b.c.'s inside the main loop for them to be

'remembered' throughout the iteration process. At the moment it

forgets them as soon as the iteration begins.

Any ideas would be most appreciated,

Geoff

Sure. At each iteration, the line u = H\u will replace the contents of u

where you first put the boundary conditions. Although

I am not sure if it is intended or not...

markus

a

N-1),

and

to

at a constant value. So that the edges of the plate are a constant

predetermined temperature.

What do you recommend i do to achieve this?

geoff

Do you really want the temperature to be 4 everywhere at time 0 but 0 at the edges? Assuming you do...

To make points stay the same temperature at all times, all you have to do is take the rows of H corresponding to those points, and

put a 1 on the diagonal and a zero everywhere else.

One more thing. The sum of every row of H should be 1. But in your case, this is only true for the rows corresponding to the

internal nodes. However, since you are specifying constant T BCs at all the edges, my point above will fix that for you.

-Paul

aul Skoczylas wrote:

a

N-1),

and

to

Ideally, my final model will enable the user to pick a point (or

points) on the grid to be a source or sink of heat. You then set an

initial temperature for the rest of the grid (an initial condition,

say the ambient temperature of the plate) and then observe the heat

transfer.

At this stage, if I can have all the boundaries at a

certain temperature that remains constant and an initial value for

the internal nodes then I'll be more than happy.

With the H matrix, wouldn't having the sum of each row of H equal to

1 mean nothing would happen during the iteration? I'm well confused

now!!

I can see what your saying about setting the H matrix so that certain

points are constant throughout the iteration, but will the

neighbouring points recognise this?

I hope I'm not being too dense about this. Feel as though I'm missing

a trick somewhere!

Thanks for your help,

Geoff

a

N-1),

and

to

Ideally, my final model will enable the user to pick a point (or

points) on the grid to be a source or sink of heat. You then set an

initial temperature for the rest of the grid (an initial condition,

say the ambient temperature of the plate) and then observe the heat

transfer.

At this stage, if I can have all the boundaries at a

certain temperature that remains constant and an initial value for

the internal nodes then I'll be more than happy.

With the H matrix, wouldn't having the sum of each row of H equal to

1 mean nothing would happen during the iteration? I'm well confused

now!!

I can see what your saying about setting the H matrix so that certain

points are constant throughout the iteration, but will the

neighbouring points recognise this?

I hope I'm not being too dense about this. Feel as though I'm missing

a trick somewhere!

Thanks for your help,

Geoff

Having the sum of each row equal to one does NOT mean nothing happens. It just means nothing wierd happens! Nothing will happen to

a specific node if the value on the diagonal is 1, and everything else is zero--but that's what you want for your constant

temperature BCs.

Exception: when you implement heat sources or sinks, I think that the rows corresponding to those nodes might not have a sum of

one.

NOTE: I don't think you're using Crank-Nicolson. That would need an extra step, which you are missing. (I've never actually done

a C-N solution of this problem--I've always used straight implicit--but I just checked a text.) This appears to be a straight

implicit method. The coefficients in the matrix don't seem quite right to me, either, but I haven't fully examined them, either. I

don't like the way that thermal diffusivity is not listed, but assumed to be (I think) 1.0. Not a very useful program if you can't

specify the actual diffusivity.

Anyway, try it! Assuming your coefficients are correct (check them again later), set the rows in your H matrix corresponding to the

edges/corners so that the diagonal is 1 and the other elements are 0. Leave the other rows alone (and those rows do already have a

sum of 1--with 1.64 on the diagonal, and four -0.16 values elsewhere). Note that 16 of your 25 rows correspond to edges--there's

only 9 real calculations going on. And those real calculations are unaffected by the trivial nature of the other rows in the

matrix.

-Paul

Use e.g. u_bc instead of u for b.c.'s.

aul Skoczylas wrote:

time 0

you

and

in

the

T

(or

set

equal

I think I've got a good idea of how to implement the b.c.'s now. The

alternative labelling will help to keep my confusion to a minimum.

The method is supposed to be crank nicholson. I borrowed the approach

from the following link:

<http://www.cs.berkeley.edu/~demmel/cs267/lecture17/lecture17.html>

not sure where i've tripped up then??!!

I haven't started looking at the coefficients yet. I've included the

diffusivity and density in my constant k. In my code it is labelled

incorrectly.

Thanks again paul and markus for all your help

time 0

you

and

in

the

T

(or

set

equal

I think I've got a good idea of how to implement the b.c.'s now. The

alternative labelling will help to keep my confusion to a minimum.

The method is supposed to be crank nicholson. I borrowed the approach

from the following link:

<http://www.cs.berkeley.edu/~demmel/cs267/lecture17/lecture17.html>

not sure where i've tripped up then??!!

I haven't started looking at the coefficients yet. I've included the

diffusivity and density in my constant k. In my code it is labelled

incorrectly.

Thanks again paul and markus for all your help

That link says (this is copied/pasted, and might not look right, depending on the font):

k k

(1 + 2*---) * U(i,j,m+1) - ----- * ( U(i-1,j,m+1) + U(i+1,j,m+1) +

h^2 2*h^2

U(i,j-1,m+1) + U(i,j+1,m+1) )

k k

= (1 - 2*---) * U(i,j,m) + ----- * ( U(i-1,j,m) + U(i+1,j,m) +

h^2 2*h^2

U(i,j-1,m) + U(i,j+1,m) )

= b(i,j,m)

This looks correct for C-N if h =dx=dy, and k=diffusivity*timestep.

The RHS of the equation is b, which is a function of U. In your code, the RHS is just U, as in your line: u = H\u (which is what

you would use for straight implicit, but for that, your coefficients wouldn't be right). You need to generate that b vector from u,

and then calculate the next time step as u=H\b. This needs to be done at every time step.

I take issue with the statement in that link that C-N is "the simplest reasonable implicit method". Straight implicit is simpler,

and at least as reasonable as the explicit method they looked at first. (More so, since like C-N, it's unconditionally stable.)

The advantage of C-N is that it's a bit more accurate, but certainly more complicated. C-N is 2nd order accurate in space and time.

The plain explicit and implicit methods are 2nd order in space, but only first order in time.

-Paul

aul Skoczylas wrote:

So if i inserted a line of code with the equation for determining b

from u, this would solve my problem? or is there more to implementing

crank nicholson than that?

Would you recommend persuing the cn method or just using straight

forward implicit?

Also if i am to use straight implicit, in what way are my

coefficients off?

Finally here is my latest code with b.c.s and i.c. sorted out. I'm

currently trying to allow for the introduction of a hot or cold spot:

%% IMPLICIT CRANK NICOLSON METHOD %%

clc;

clear all;

% define the constants for the problem

M = 25; % number of time steps

L = 1; % length and width of plate

k = 0.01; % constant

N = 30 % number of elements in plate (NxN)

h = L/(N-1); % grid spacing

x = [0:h:L]; % vector of x values

y = [0:h:L]; % vector of y values

%

[xx yy] = meshgrid(x,y);

%

% set constants a and c

a = (1 + (2*k)/(h^2));

c = (-k/(2*(h^2)));

%

% set the matrix for H

A = a*eye(N) + c*diag(ones(1, N-1), 1) + c*diag(ones(1, N-1),

-1);

B = c*eye(N);

C = zeros(N);

% Create the matrix H, composed of sub-matrices A,B and C

Ap = eye(N);

Bp = diag(ones(1, N-1), 1) + diag(ones(1, N-1), -1);

Cp = ones(N) - Ap - Bp;

%

% Define the H matrix using KRON to "replace" 1's in Ap, Bp, and Cp

with copies of A, B, and C

H = kron(Ap, A) + kron(Bp, B) + kron(Cp, C);

%

% % set up initial conditions for plate

a = 4; % the initial condition value for the plate

u = a*ones(N,N);

%

% set the boundary conditions for the plate

top = 4;

u(1,:) = top; % boundary condition for top of plate

left = 4;

u(:,1) = left; % boundary condition for left of plate

bottom = 4;

u(N,:) = bottom;% boundary condition for bottom

right = ;

u(:,N) = right; % boundary condition for right of plate

%

% set hot or cold spot

source_sink = 2;

k = 5; l = 10;

x_posn = k:l;

p = 5; q = 10;

y_posn = p:q;

u(y_posn,x_posn) = source_sink;

%

u = u(:);

%

for m = 1:M

u = H\u;

% replace the b.c. values

u(1:N) = left; % left b.c.

u(((2*N):N:N*N)) = bottom; % bottom b.c.

u((N+1):N:(N*N)) = top; % top b.c.

u((((N-1)*N)+1):(N^2)) = right; % right b.c.

% replace source/sink value

u((N*(k-1)+p):(N*(k-1)+q)) = source_sink;

u((N*(k)+p):(N*(k)+q)) = source_sink;

u((N*(k+1)+p):(N*(k+1)+q)) = source_sink;

u((N*(k+2)+p):(N*(k+2)+q)) = source_sink;

u((N*(k+3)+p):(N*(k+3)+q)) = source_sink;

u((N*(l-1)+p):(N*(l-1)+q)) = source_sink;

end

% create a surf plot of the temperature (z-axis) against grid

position

axis manual; axis([0 L 0 L 0 a]),% pause;

drawnow

uu = reshape(u,N,N)

surf(xx,yy,uu), colorbar

What do you think? The part for replacing source/sink values is a

little crude. I couldn't figure out how to define the set temperature

for each column of u and therefore just did each in turn.

Thanks again,

Geoff

So if i inserted a line of code with the equation for determining b

from u, this would solve my problem? or is there more to implementing

crank nicholson than that?

Would you recommend persuing the cn method or just using straight

forward implicit?

Also if i am to use straight implicit, in what way are my

coefficients off?

Finally here is my latest code with b.c.s and i.c. sorted out. I'm

currently trying to allow for the introduction of a hot or cold spot:

%% IMPLICIT CRANK NICOLSON METHOD %%

clc;

clear all;

% define the constants for the problem

M = 25; % number of time steps

L = 1; % length and width of plate

k = 0.01; % constant

N = 30 % number of elements in plate (NxN)

h = L/(N-1); % grid spacing

x = [0:h:L]; % vector of x values

y = [0:h:L]; % vector of y values

%

[xx yy] = meshgrid(x,y);

%

% set constants a and c

a = (1 + (2*k)/(h^2));

c = (-k/(2*(h^2)));

%

% set the matrix for H

A = a*eye(N) + c*diag(ones(1, N-1), 1) + c*diag(ones(1, N-1),

-1);

B = c*eye(N);

C = zeros(N);

% Create the matrix H, composed of sub-matrices A,B and C

Ap = eye(N);

Bp = diag(ones(1, N-1), 1) + diag(ones(1, N-1), -1);

Cp = ones(N) - Ap - Bp;

%

% Define the H matrix using KRON to "replace" 1's in Ap, Bp, and Cp

with copies of A, B, and C

H = kron(Ap, A) + kron(Bp, B) + kron(Cp, C);

%

% % set up initial conditions for plate

a = 4; % the initial condition value for the plate

u = a*ones(N,N);

%

% set the boundary conditions for the plate

top = 4;

u(1,:) = top; % boundary condition for top of plate

left = 4;

u(:,1) = left; % boundary condition for left of plate

bottom = 4;

u(N,:) = bottom;% boundary condition for bottom

right = ;

u(:,N) = right; % boundary condition for right of plate

%

% set hot or cold spot

source_sink = 2;

k = 5; l = 10;

x_posn = k:l;

p = 5; q = 10;

y_posn = p:q;

u(y_posn,x_posn) = source_sink;

%

u = u(:);

%

for m = 1:M

u = H\u;

% replace the b.c. values

u(1:N) = left; % left b.c.

u(((2*N):N:N*N)) = bottom; % bottom b.c.

u((N+1):N:(N*N)) = top; % top b.c.

u((((N-1)*N)+1):(N^2)) = right; % right b.c.

% replace source/sink value

u((N*(k-1)+p):(N*(k-1)+q)) = source_sink;

u((N*(k)+p):(N*(k)+q)) = source_sink;

u((N*(k+1)+p):(N*(k+1)+q)) = source_sink;

u((N*(k+2)+p):(N*(k+2)+q)) = source_sink;

u((N*(k+3)+p):(N*(k+3)+q)) = source_sink;

u((N*(l-1)+p):(N*(l-1)+q)) = source_sink;

end

% create a surf plot of the temperature (z-axis) against grid

position

axis manual; axis([0 L 0 L 0 a]),% pause;

drawnow

uu = reshape(u,N,N)

surf(xx,yy,uu), colorbar

What do you think? The part for replacing source/sink values is a

little crude. I couldn't figure out how to define the set temperature

for each column of u and therefore just did each in turn.

Thanks again,

Geoff

Geoff Halstead" < XXXX@XXXXX.COM > wrote in message news: XXXX@XXXXX.COM ...

I think that would do it.

It depends what your goal is. If you want to do it as simply as possible, I'd forget C-N, and go straight implicit. But if you're

going to be making a code which is going to be running a huge simulation and computer time and accuracy are both very important,

than C-N is better. (But if computer time is very important, there's a whole lot you need to do, because this is currently very

inefficient.)

Where you now have (again, assuming k=diffusivity*timestep and h=dx=dy) for C-N

a = (1 + (2*k)/(h^2));

c = (-k/(2*(h^2)));

You would need for fully implicit

a = (1 + (4*k)/(h^2));

c = (-k/(h^2));

You really should try to understand where these come from! They are from a simple discretization of the heat conduction equation.

You simply insert numerical derivatives into the equation to replace the derivatives that are there. The difference between

explicit, implicit and C-N is at what time you evaluate the spatial derivatives.

Definitely crude! I really don't like the idea of replacing the temperatures at every step like that. But I've already made my

recommendation for how to do this instead... And since your sink/source is really only another constant temperature BC, the same

would work for it. Why make a change at every time step if you only need to do it once? Very inefficient!

I also get very concerned when you reuse variables. In your original code, you reused a, and that caused me some confusion. In

this, you still have that, but now you're also reusing k, and that is very dangerous--especially if you decide to continue with C-N,

since you're going to need the original k value!

Note that the improper implementation of C-N is causing errors. It looks to me like this plate is at 4 everywhere at the start.

The edges are held at 4 throughout, while a block in the middle is held at 2 after time 0. Therefore there should never be any

temperature on the plate which is not between 2 and 4. But when I run your code, the temperature is less than 2 everywhere except

in the regions where it is specified as 2 or 4.

Aside from all that, here's a way to speed it up:

The line u=H\u can be replaced with:

[LL,UU,PP]=lu(H);

u=UU\(LL\(PP*u));

(Of course, if you fix the C-N implementation, the u on the right will have to change to b.)

The trick here is that since H never changes, the first line can be put outside the loop, and only the second part goes inside the

loop. On my computer, doing this dropped the execution time by a factor of 10!

I think that would do it.

It depends what your goal is. If you want to do it as simply as possible, I'd forget C-N, and go straight implicit. But if you're

going to be making a code which is going to be running a huge simulation and computer time and accuracy are both very important,

than C-N is better. (But if computer time is very important, there's a whole lot you need to do, because this is currently very

inefficient.)

Where you now have (again, assuming k=diffusivity*timestep and h=dx=dy) for C-N

a = (1 + (2*k)/(h^2));

c = (-k/(2*(h^2)));

You would need for fully implicit

a = (1 + (4*k)/(h^2));

c = (-k/(h^2));

You really should try to understand where these come from! They are from a simple discretization of the heat conduction equation.

You simply insert numerical derivatives into the equation to replace the derivatives that are there. The difference between

explicit, implicit and C-N is at what time you evaluate the spatial derivatives.

Definitely crude! I really don't like the idea of replacing the temperatures at every step like that. But I've already made my

recommendation for how to do this instead... And since your sink/source is really only another constant temperature BC, the same

would work for it. Why make a change at every time step if you only need to do it once? Very inefficient!

I also get very concerned when you reuse variables. In your original code, you reused a, and that caused me some confusion. In

this, you still have that, but now you're also reusing k, and that is very dangerous--especially if you decide to continue with C-N,

since you're going to need the original k value!

Note that the improper implementation of C-N is causing errors. It looks to me like this plate is at 4 everywhere at the start.

The edges are held at 4 throughout, while a block in the middle is held at 2 after time 0. Therefore there should never be any

temperature on the plate which is not between 2 and 4. But when I run your code, the temperature is less than 2 everywhere except

in the regions where it is specified as 2 or 4.

Aside from all that, here's a way to speed it up:

The line u=H\u can be replaced with:

[LL,UU,PP]=lu(H);

u=UU\(LL\(PP*u));

(Of course, if you fix the C-N implementation, the u on the right will have to change to b.)

The trick here is that since H never changes, the first line can be put outside the loop, and only the second part goes inside the

loop. On my computer, doing this dropped the execution time by a factor of 10!

aul Skoczylas wrote:

determining

straight

N-1),

and

a

repeated ones. The new line for u = H\u really helps!

I'm now trying to change my code from crank nicholson to a simple

implicit scheme by altering the H matrix.

I'm going to impliment you idea for b.c.s later on, when I get the

code working and doing what I want it to.

cheers,

Geoff

determining

straight

N-1),

and

a

repeated ones. The new line for u = H\u really helps!

I'm now trying to change my code from crank nicholson to a simple

implicit scheme by altering the H matrix.

I'm going to impliment you idea for b.c.s later on, when I get the

code working and doing what I want it to.

cheers,

Geoff

aul Skoczylas wrote:

determining

straight

N-1),

and

a

I was just wondering which implicit scheme you would recommend for me

to use. I eventually plan to translate this problem into 3D, so I'm

not sure whether this will have a bearing on the decision.

How easy is it to include the b(m) for c-nicholson in the code i have

at the moment?

cheers, geoff

determining

straight

N-1),

and

a

I was just wondering which implicit scheme you would recommend for me

to use. I eventually plan to translate this problem into 3D, so I'm

not sure whether this will have a bearing on the decision.

How easy is it to include the b(m) for c-nicholson in the code i have

at the moment?

cheers, geoff

gt;> Definitely crude! I really don't like the idea of replacing

the

C-N,

when

right

can

a

just a quick one. I've altered my code so it is now crank nicholson

implicit. I'm now having a few problems with my constant k. When its

very small the results are what i expect but if i set k to anything

above 0.00001 it starts doing strange things on the corners of the

boundaries. what do you think is going wrong? here's the latest:

% define the constants for the problem

M = 100; % number of time steps

L = 1; % length and width of plate

k = 5; % constant

N = 30; % number of elements in plate (NxN)

h = L/(N-1); % grid spacing

x = [0:h:L]; % vector of x values

y = [0:h:L]; % vector of y values

%

[xx yy] = meshgrid(x,y);

%

% set constants a and c

a = (1 + (2*k)/(h^2));

c = (-k/(2*(h^2)));

d = (1 - (2*k)/(h^2));

e = (k/(2*(h^2)));

%

% set the matrix for H

A = a*eye(N) + c*diag(ones(1, N-1), 1) + c*diag(ones(1, N-1),

-1);

B = c*eye(N);

D = d*eye(N) + e*diag(ones(1, N-1), 1) + e*diag(ones(1, N-1),

-1);

E = e*eye(N);

C = zeros(N);

% Create the matrix H, composed of sub-matrices A,B and C

Ap = eye(N);

Bp = diag(ones(1, N-1), 1) + diag(ones(1, N-1), -1);

Dp = eye(N);

Ep = diag(ones(1, N-1), 1) + diag(ones(1, N-1), -1);

Cp = ones(N) - Ap - Bp;

%

% Define the H_1 matrix using KRON to "replace" 1's in Ap, Bp, and Cp

with copies of A, B, and C

H_1 = kron(Ap, A) + kron(Bp, B) + kron(Cp, C);

% Define the H_2 matrix using KRON to "replace" 1's in Ap, Bp, and Cp

with copies of A, B, and C

H_2 = kron(Dp, D) + kron(Ep, E) + kron(Cp, C);

%

% set up initial conditions for plate

i_c = 4; % the initial cond value for the plate

u = i_c*ones(N,N);

%

% set the boundary conditions for the plate

top = 4;

u(1,:) = top; % bc for top of plate

left = 4;

u(:,1) = left; % bc for left side of plate

bottom = 4;

u(N,:) = bottom; % bc for bottom of plate

right = 4;

u(:,N) = right; % bc for right side of plate

%

% set hot or cold spot

source_sink = 2;

r = 5; s = 10;

x_posn = r:s;

p = 5; q = 10;

y_posn = p:q;

u(y_posn,x_posn) = source_sink;

%

u = u(:);

%

[SS_1,UU_1,PP_1] = lu(H_1);

%

for m = 1:M

b = H_2*u;

u = UU_1\(SS_1\(PP_1*b));

% replace the b.c. values

u(1:N) = left; % left b.c.

u(((2*N):N:N*N)) = bottom; % bottom b.c.

u((N+1):N:(N*N)) = top; % top b.c.

u((((N-1)*N)+1):(N^2)) = right; % right b.c.

% replace source/sink value

u((N*(r-1)+p):(N*(r-1)+q)) = source_sink;

u((N*(r)+p):(N*(r)+q)) = source_sink;

u((N*(r+1)+p):(N*(r+1)+q)) = source_sink;

u((N*(r+2)+p):(N*(r+2)+q)) = source_sink;

u((N*(r+3)+p):(N*(r+3)+q)) = source_sink;

u((N*(s-1)+p):(N*(s-1)+q)) = source_sink;

end

geoff

the

C-N,

when

right

can

a

just a quick one. I've altered my code so it is now crank nicholson

implicit. I'm now having a few problems with my constant k. When its

very small the results are what i expect but if i set k to anything

above 0.00001 it starts doing strange things on the corners of the

boundaries. what do you think is going wrong? here's the latest:

% define the constants for the problem

M = 100; % number of time steps

L = 1; % length and width of plate

k = 5; % constant

N = 30; % number of elements in plate (NxN)

h = L/(N-1); % grid spacing

x = [0:h:L]; % vector of x values

y = [0:h:L]; % vector of y values

%

[xx yy] = meshgrid(x,y);

%

% set constants a and c

a = (1 + (2*k)/(h^2));

c = (-k/(2*(h^2)));

d = (1 - (2*k)/(h^2));

e = (k/(2*(h^2)));

%

% set the matrix for H

A = a*eye(N) + c*diag(ones(1, N-1), 1) + c*diag(ones(1, N-1),

-1);

B = c*eye(N);

D = d*eye(N) + e*diag(ones(1, N-1), 1) + e*diag(ones(1, N-1),

-1);

E = e*eye(N);

C = zeros(N);

% Create the matrix H, composed of sub-matrices A,B and C

Ap = eye(N);

Bp = diag(ones(1, N-1), 1) + diag(ones(1, N-1), -1);

Dp = eye(N);

Ep = diag(ones(1, N-1), 1) + diag(ones(1, N-1), -1);

Cp = ones(N) - Ap - Bp;

%

% Define the H_1 matrix using KRON to "replace" 1's in Ap, Bp, and Cp

with copies of A, B, and C

H_1 = kron(Ap, A) + kron(Bp, B) + kron(Cp, C);

% Define the H_2 matrix using KRON to "replace" 1's in Ap, Bp, and Cp

with copies of A, B, and C

H_2 = kron(Dp, D) + kron(Ep, E) + kron(Cp, C);

%

% set up initial conditions for plate

i_c = 4; % the initial cond value for the plate

u = i_c*ones(N,N);

%

% set the boundary conditions for the plate

top = 4;

u(1,:) = top; % bc for top of plate

left = 4;

u(:,1) = left; % bc for left side of plate

bottom = 4;

u(N,:) = bottom; % bc for bottom of plate

right = 4;

u(:,N) = right; % bc for right side of plate

%

% set hot or cold spot

source_sink = 2;

r = 5; s = 10;

x_posn = r:s;

p = 5; q = 10;

y_posn = p:q;

u(y_posn,x_posn) = source_sink;

%

u = u(:);

%

[SS_1,UU_1,PP_1] = lu(H_1);

%

for m = 1:M

b = H_2*u;

u = UU_1\(SS_1\(PP_1*b));

% replace the b.c. values

u(1:N) = left; % left b.c.

u(((2*N):N:N*N)) = bottom; % bottom b.c.

u((N+1):N:(N*N)) = top; % top b.c.

u((((N-1)*N)+1):(N^2)) = right; % right b.c.

% replace source/sink value

u((N*(r-1)+p):(N*(r-1)+q)) = source_sink;

u((N*(r)+p):(N*(r)+q)) = source_sink;

u((N*(r+1)+p):(N*(r+1)+q)) = source_sink;

u((N*(r+2)+p):(N*(r+2)+q)) = source_sink;

u((N*(r+3)+p):(N*(r+3)+q)) = source_sink;

u((N*(s-1)+p):(N*(s-1)+q)) = source_sink;

end

geoff

Without looking at the code, here's a thought:

Your k constant is really thermal diffusivity multiplied by the time step. I don't know what material you're trying to model, but

the diffusivity of steel is on the order of 0.00002 ms or so, depending on the grade. Diffusivity of wood is 0.0000001 ms or

so.

So looking at steel, if your k is 0.00001, then your time step is 0.5 seconds. That's a pretty big time step. In wood, it's much

larger still!

Both C-N and straight implicit are unconditionally stable at all time step sizes--they don't blow up if your time step is too

large--but they aren't unconditionally accurate. You will get better accuracy with small time steps. You are just seeing

inaccuracies related to large time steps. These inaccuracies should settle out with time--you should approach the correct steady

state solution eventually. (But like I said, I didn't look at your code--there certainly could be other errors in there!)

I'd recommend making a simplified problem to which you know what the transient solution should be, and making sure your model is

correct with it, before you trust it on a complex problem. (Unfortunately, transient problems of this nature usually require a sum

of an infinite series to calculate analytically; often each term in the infinite series requires a solution to a transcendental

equation.)

-Paul

Similar Threads

1. Heat equation with neumann boundary conditions

I am looking to adapt a simple 1 dimensional heat quation code. The code I am currently using makes use of dirichlet boundary conditions, that is to say the end temperatures of the rod are set at a specific temperature. I would however like to use neumann b.c.'s which effectively model heat flow in or out of the end of the rod. I have looked for a method of applying a b.c. which is time dependent but cannotfigure out how to apply it to my code. Here is my matlab code, which uses the 2nd order implicit crank nicholson method: %% IMPLICIT CRANK NICOLSON METHOD %% clc; clear all; % define the constants for the problem L = 1; % length of bar k = 0.5; % thermal conductivity of bar N = 3; % number of elements in bar J = 200; % number of time steps dx = L/(N+1); x = [0:dx:L]; % vector of x values s = 0.6*k; % s is the step size % set the matrix for u and the initial conditions Id = eye(N+2); % sets up a zeros matrix of the desired dimensions S = ones(N+2); % S = tril(S,1); S = triu(S,-1); S = S - 3*Id; S(1,1) = 0; S(1,2) = 0; S(N+2,N+1) = 0; S(N+2,N+2) = 0 ; A = Id +s*S; A = inv(Id - s*S)*A; % Crank nihcolson part C =A' % % % set up initial conditions for bar u = ones(1,N+2); u = u.*0; %input('initial condition = '); % set the boundary conditions for the bar u(1) = 1; %input('left boundary condition = '); % boundary condition for left hand end of bar u(N+2) = 0 %input('right boundary condition = '); % boundary condition for right hand end of bar for j = 1:J u = u*C; drawnow; plot(x,u); end Any ideas as to how to create time dependent b.c.s would be most appreciated, Geoff

2. how to specify boundary condition for heat equation - MATLAB

3. specify boundary condition for heat equation

How to specific boundary condition for heat equation that at one side the B.C. takes the form of dU/dt=0? thanks!

4. Boundary conditions definition for heat equation PDE Toolbox - MATLAB

5. 2nd Order wave equation in 2D, SPECTRAL METHODS, implementation of BOUNDARY CONDITIONS

The following script (from Trefethen, Spectral Methods in MATLAB) solves the 2nd order wave equation in 2 dimensions (u_tt=u_xx+u_yy) using spectral methods, Fourier for x and Chebyshev for y direction. On its rectangular domain, the equation is subject to Neumann boundary conditions along the sides ( u_y(x,+-1,t) =0 ) and periodic boundary conditions at the ends ( u (+3,y,t)= u(-3,y,t) ). The calculation in BC = -Dy([1 Ny+1],[1 Ny+1])\Dy([1 Ny+1],2:Ny) is a (pre)solution of a system of equations (I suppose the one coming up in vv([1 Ny+1],:) = BC*vv(2:Ny,:) near the end) and Dy's are Chebyshev spectral differentiation matrices in y direction. This BC is then used near the end of the code in vv([1 Ny+1],:) = BC*vv(2:Ny,:) where the Boundary conditions are actually implemented. The first and last rows of the solution vv (corresponding to y=+-1) are being updated and the Neumann conditions are being used. I don't seem to be able to see what is EXACTLY happending there. How is vv([1 Ny+1],:) = BC*vv(2:Ny,:) related to u_y(x,+- 1,t)=0? Thanks, Omid. PS. This script requires the function cheb() which I'm also posting, in case anyone wanted to actually try the code. % p37.m - 2D "wave tank" with Neumann BCs for |y|=1 % x variable in [-A,A], Fourier: A = 3; Nx = 50; dx = 2*A/Nx; x = -A+dx*(1:Nx)'; D2x = (pi/A)^2*toeplitz([-1/(3*(dx/A)^2)-1/6 ... .5*(-1).^(2:Nx)./sin((pi*dx/A)*(1:Nx-1)/2).^2]); % y variable in [-1,1], Chebyshev: Ny = 15; [Dy,y] = cheb(Ny); D2y = Dy^2; BC = -Dy([1 Ny+1],[1 Ny+1])\Dy([1 Ny+1],2:Ny); % Grid and initial data: [xx,yy] = meshgrid(x,y); vv = exp(-8*((xx+1.5).^2+yy.^2)); dt = 5/(Nx+Ny^2); vvold = exp(-8*((xx+dt+1.5).^2+yy.^2)); % Time-stepping by leap frog formula: clf, plotgap = round(2/dt); dt = 2/plotgap; for n = 0:2*plotgap t = n*dt; if rem(n+.5,plotgap)<1 subplot(3,1,n/plotgap+1), mesh(xx,yy,vv), view(-10,60) axis([-A A -1 1 -0.15 1]), colormap(1e-6*[1 1 1]); text(-2.5,1,.5,['t = ' num2str(t)],'fontsize',18) set(gca,'ztick',[]), grid off, drawnow end vvnew = 2*vv - vvold + dt^2*(vv*D2x +D2y*vv); vvold = vv; vv = vvnew; vv([1 Ny+1],:) = BC*vv(2:Ny,:); end % CHEB compute D = differentiation matrix, x = Chebyshev grid function [D,x] = cheb(N) if N==0, D=0; x=1; return, end x = cos(pi*(0:N)/N)'; c = [2; ones(N-1,1); 2].*(-1).^(0:N)'; X = repmat(x,1,N+1); dX = X-X'; D = (c*(1./c)')./(dX+(eye(N+1))); % off-diagonal entries D = D - diag(sum(D')); % diagonal entries

6. 1 D Heat equation solving by Crank Nicolson method - MATLAB

Hi guys I need a bit from help, have to do a program using Crank Nicholson the equation of the heat for a plate. The program serious seemed to this one, that it is of a bar. Can someone orientate me on since doing it? Thank you function MWD=calorcrank(m,N,L,c,T) h=L/m;k=T/N; lan=c^2*k/h^2; %red espacial for i=1:m-1 x(i)=i*h; end %matriz A (coeficientes del sistema que tenemos que resolver) %es una matriz ftr tridiagonal A=zeros(m-1); A(1,1)=2+2*lan; A(1,2)=-lan; A(m-1,m-2)=-lan; A(m-1,m-1)=2+2*lan; for i=2:m-2 A(i,i-1)=-lan; A(i,i)=2+2*lan; A(i,i+1)=-lan; end %matriz B (coeficientes del sistema que tenemos que resolver) %es una matriz ftr tridiagonal B=zeros(m-1); B(1,1)=2-2*lan; B(1,2)=lan; B(m-1,m-2)=lan; B(m-1,m-1)=2-2*lan; for i=2:m-2 B(i,i-1)=lan; B(i,i)=2-2*lan; B(i,i+1)=lan; end %Calculo de W0, le llamamos W %equivalentemente W=ff(x) for i=1:m-1 W(i)=ff(x(i)); end W=W'; cont=1; MW(cont,:)=W'; %Esquema iterativo for j=1:N %calculo del vector variable en j, bb bb(1)=k*F(x(1),(j-1/2)*k)+lan*a(j*k)+lan*a((j-1)*k); bb(m-1)=k*F(x(m-1),(j-1/2)*k)+lan*b(j*k)+lan*b((j-1)*k); for i=2:m-2 bb(i)=k*F(x(i),(j-1)*k); end cc=B*W+bb'; %terminos independientes %solucion del sistema de ecuaciones lineales AW=cc %en este caso lo resolvemos por el metodo de eliminacion %gaussiana que es el que tiene MATLAB por defecto W=A\cc W=A\cc; cont=cont+1; MW(cont,:)=W'; end t=0:k:T; vi=a(t); vf=b(t); MWD(:,1)=vi'; MWD(:,m+1)=vf'; MWD(:,2:m)=MW; En otros ficheros Mfiles: >> edit a.m function y=a(t) y=exp(-t).*(1-t); >> edit b.m function y=b(t) y=0; >> edit F.m function y=F(x,t) y=exp(-t).*(2-t).*(x-1); >> edit ff.m function y=ff(x) y=1-x+sin(pi.*x);