Find temperature vector using finite difference method

I am aware that there is a very similar question posted here that asks for the same thing, however I don't really understand what they have done nor does it seem applicable to my problem. Just like the other problem I have the differential equation:

where T0=300K, TL = 400K and k=2.5.

And I am given that the second derivative with central difference is:

Further information given is that N = 250. I had previously determined that matrix A would be: [2 -1 0; -1 2 -1; 0 -1 2] when N=4. However, I am unsure if this is the matrix that I should use or if I should use a general matrix such as:

diagonal = [1/h^2; 2/h^2*ones(N-1,1); 1/h^2];
A = diag( diagonal)

for i = 2:N+1
    A(i-1,i) = -1/h^2; A(i,i-1)=-1/h^2;
end

A(1,2)=0; A(N+1,N)=0;

If I use the general matrix and the following code:

L = 2;
N = 250;
T0 = 300; % boundary condition
Tend = 400; % boundary condition
h = L/N;
xj = 0:h:1;
Q=@(x)(300*exp(-(xj-L/2)).^2);

diagonal = [1/h^2; 2/h^2*ones(N-1,1); 1/h^2];
A = diag( diagonal);
for i = 2:N+1
    A(i-1,i)=1/h^2; A(i,i-1)=1/h^2;
end
A(1,2)=0; A(N+1,N)=0;

b = [T0/h^2; Q; Tend/h^2];

T=A\b
plot(xj,T)

I get an error message at Q saying I have too many input arguments. I don't know if that is the correct way to find Q for the different x values and then how to code that be should contain all these values. Could anyone please tell me how to define Q and if the rest of the code (using the general matrix etc.) is correct. Thank you!

Answer 1

[Alan Stevens]

Try this

L = 2;
N = 250;
T0 = 300; % boundary condition
Tend = 400; % boundary condition
h = L/N;
xj = 0:h:L;  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% L not 1
Q=@(x)300*exp(-(x-L/2).^2);  %%%%%%%%%%%%%%%%%  x not xj.  Also you had wrong bracket positions

diagonal = [1/h; 2/h^2*ones(N-1,1); 1/h];
A = diag( diagonal);
for i = 2:N+1
    A(i-1,i)=-1/h^2; A(i,i-1)=-1/h^2;  %%%%%%%%%%%%%%  signs
end
A(1,2)=0; A(N+1,N)=0;

b = [T0/h; Q(xj(2:end-1))'; Tend/h];  %%%%%%%%%%%%%%% Q(xj(2_end-1))' not just Q

T=A\b;
plot(xj,T),grid

to get

[Alan Stevens]

Also, you don't seem to have included a thermal conductivity (k). (Unless you deliberately set it to 1, of course).

[Katara So]

Thank you!

I didn't think that k mattered for this problem but it is given that k = 2.5. Does that change the code?

[Katara So]

Also, may I ask why Q becomes

Q(xj(2:end-1))'

[Alan Stevens]

1. k will only affect the magnitude of the response.

2. The transpose is there to make the whole thing a column vector. The 2:end-1 is to make b the right size. The first and last temperatures are fixed, so you don't need Q for them (an alternative would be to only have calculated Q for a reduced number of nodes in the first place).

[Katara So]

Okay, so I don't need to include k here?

[Alan Stevens]

You don't need k to make the routine work. However, if the other parameters have realistic values it seems odd to leave out k.

[Katara So]

In the problem I linked in my original question they had that k = 2 + x/7 which was then used in their code in some way I didn't understand. So if I want to add that k = 2.5 or if k is some equation as above, should I change my code to look like the one I linked? Or is there some other way to use k? Truly appreciate all the help!

[Alan Stevens]

The following shows one way of incorporating k

L = 2;
N = 250;
T0 = 300; % boundary condition
Tend = 400; % boundary condition
h = L/N;
xj = 0:h:L;  
Q=@(x)300*exp(-(x-L/2).^2);  
%k = 2.5;              %%%%%%%% Select this or the next line as desired
k = 2+xj(2:end-1)'/7;  %%%%%%%%
diagonal = [1/h; 2/h^2*ones(N-1,1); 1/h];
A = diag( diagonal);
for i = 2:N+1
    A(i-1,i)=-1/h^2; A(i,i-1)=-1/h^2;  
end
A(1,2)=0; A(N+1,N)=0;

b = [T0/h; Q(xj(2:end-1))'./k; Tend/h]; %%%%%%%  Notice Q./k now included.

T=A\b;
plot(xj,T),grid

[Katara So]

Ohh okay! Thank you very much!

[Alan Stevens]

On second thoughts the situation is a little more complicated if k is a function of x. In that case we have to go back to the original equations which are now

though the difference will be small if the change in k is small over the range of x.

[Saqib Zafar]

can you answer the task 2. plz help me

[Alan Stevens]

@Saqib Zafar. ??? I've no idea what task 2 is!!!!

[Saqib Zafar]

[Saqib Zafar]

[Saqib Zafar]

plz help me in this task..plz