Find temperature vector using finite difference method

4 visualizaciones (últimos 30 días)
Katara So
Katara So el 18 de Sept. de 2020
Comentada: Saqib Zafar el 28 de Sept. de 2020
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!

Respuesta aceptada

Alan Stevens
Alan Stevens el 18 de Sept. de 2020
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
  15 comentarios
Saqib Zafar
Saqib Zafar el 28 de Sept. de 2020
plz help me in this task..plz

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Programming en Help Center y File Exchange.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by