# Where is the error on my algorithm?

1. Jan 2, 2012

### loukoumas

Hello everyone, happy new year!

I have already got much help from your forum on solving with numerical methods heat diffusion problems. I am afraid a little that in this thread i ask too much... But maybe you can still help me.

I use "Crank-Nicolson" method to calculate the temperatures,as they change on time, on 1m length aluminum rod (which is assumed to be culindric insolated , so diffusion acts on one dimension only). According to this method, i use 11 points to divide my rod in 10 intervals and i use the equation:
T(i point)=a*E(i point)+b*[E(i-1 point)+E(i+1 point)+T(i-1 point)+T(i+1 point)

E(i) is the temperature of "i" point but at the previous point of time (i have divided time in 1 sec intervals so if E(i) is the temperature on 0 seconds then T(i) is on 1 seconds ). a,b are known constants and also T(1)=E(1)=0, T(11)=E(11)=100.
So if i know T(0m,t) T(1m,t) and T(x,0 s) , i wright my equation for all my points, except i=1 and i=11 where temperatures are known (i could wright i=0 and i=10 but i use "i" to create a 1-column matrix in matlab so i can't use i=0).

T(2)=A*E(2)+B*[E(1)+E(3)+T(1)+T(3)]
T(3)=A*E(3)+B*[E(2)+E(4)+T(2)+T(4)]
T(4)=A*E(4)+B*[E(3)+E(5)+T(3)+T(5)]
:
:
T(10)=A*E(10)+B*[E(9)+E(11)+T(9)+T(11)]

this is a 9X9 system and i follow two ways to solve it
1. I take all the unknown T(2-10) on the left part and the known E and T(1),T(11) on the right part. Then we use E and T as 1-column matrixes that are multiplyed with 9X9 matrixes. We conclude in something like this: [M]*[T]=a*[E]+[N]*[E]+b*[C] where T,E,C are 9X1 and M,N are 9X9 matixes. We multiply with the inverted M matrix and we get T

2.Back to the system of equations again. If we knew T(3) we would calculate T(2). Having already calculated T(2) we could calculate T(3) if we knew T(4). Here is what we do in this method. We give a random value at T(3),T(4)...T(10) on the right part,so we get a (incorrect)value for T(2),T(3)...T(10) on the left part. If we repeat this calculations we see that T(2),T(3)...T(10) on the left part are converging to a value, which is their right value.

When T(i) are calculated for the 1st second, we use them in E(i)'s place so we calculate T(i) for the 2nd second. So if we want to calculate temperatures for the 50th second we repeat the calculation with new every time E(i) for 50 times

I have made algorithms in matlab for both ways so i calculate T(i) ( but i send just one)

Where is my problem??
if T(0,t)=0 ⁰C T(L,t)=100 ⁰C and T(x,0)=20 ⁰C both algorithms give the same results

BUT
Trying to solve a little difrent problem with the right end of the rod insolated (so T(0,t)=0 ⁰C,T(x,0)=20 ⁰C but T(L,t)≠100 ⁰C), which means according to "Crank-Nicolson"

T(11)=(4/3)*T(10) - (1/3)*T(9)

T(11) is not constant now!

the second way doesn't give correct results. In fact both algorithms give almost the same results for the first 1000 seconds and then the differences are getting biger and biger specially for the last points of the rod(T(11),T(10),T(9))

After the necessary seconds are past, the rod must be on 0 ⁰C
But only the algorithm that uses M,N and C matrixes gives correct results
I need to find the problem on the convergent values method for the insolated right end problem

I send you one algorithm (the one for the convergent values method) for the T(0,t)=0 ⁰C T(L,t)=100 ⁰C and T(x,0)=20 ⁰C problem
and two algorithms for the T(0,t)=0 ⁰C, T(x,0)=20 ⁰C and T(11)=(4/3)*T(10) - (1/3)*T(9) which is the insolated end problem (the "matrixes 2" algorithm gives correct results). They are simple enough to follow them i think!
( i tried to send four files, two algorithms for each problem, but when i compressed them the file was invalid on the attachements. So i chose not to send the method with the matrixes for the easy,"constant temperatures on the ends" problem )

I think it is too much what i m asking for, and maybe it's not a very intresting problem to spend your time on! But if you think you can help i would be very lucky!( in a two dimensions heat diffusion problem, convergent values method is much more useful,that's why i want to find where the problem is)

Thank you very much!!