Finite difference method problem with solving an equation

Trying to use Finite difference method, to write the equation in AT = b matrices. But I don't know how to write FDM on that type of equation, please see image.

Answer 1

[Torsten]

(k(x_(i+1/2))*(T(x_(i+1))-T(x_(i)))/(x_(i+1)-x_(i))-k(x_(i-1/2))*(T(x_(i))-T(x_(i-1)))/(x_(i)-x_(i-1)))/(x_(i+1/2)-x_(i-1/2)) = Q(x_i)

with

x_(i+1/2)=(x_(i+1)+x_(i))/2
x_(i-1/2)=(x_(i)+x_(i-1))/2

Best wishes

Torsten.

Answer 2

[Torsten]

[d/dx(k*dT/dx)](x) ~ (k(x+h/2)*(dT/dx)(x+h/2) - k(x-h/2)*(dT/dx)(x-h/2))/h

Approximations for (dT/dx)(x+h/2) and dT/dx(x-h/2) are

(dT/dx)(x+h/2) ~ (T(x+h)-T(x))/h 
(dT/dx)(x-h/2) ~ (T(x)-T(x-h))/h

Now insert in the first equation.

The finite-difference approximation in my first response was more general because it took into account non-equidistant grids (i.e. h is not fixed over the complete interval). But note that I missed the minus-sign in front of the approximaton for d/dx(k*dT/dx).

Best wishes

Torsten.

Answer 3

[Torsten]

-(k(x_(i+1/2))*y_(i+1)-(k(x_(i+1/2))+k(x_(i-1/2)))*y_(i)+k(x_(i-1/2))*y_(i-1))/h^2 = Q(x_(i))

with

x_(i+1/2)=(x_(i)+x_(i+1))/2
x_(i-1/2)=(x_(i-1)+x_(i))/2

for i=2,...,(n-1)

y(1) and y(n) are given by your boundary conditions.

So to incorporate the boundary conditions, you can add the equations

y(1) = 300
y(n) = 410

to the system from above.

The complete system to be solved then looks like

y(1) = 300
-(k(x_(i+1/2))*y_(i+1)-(k(x_(i+1/2))+k(x_(i-1/2)))*y_(i)+k(x_(i-1/2))*y_(i-1))/h^2 = Q(x_(i))  (i=2,...,n-1)
y(n) = 410

Best wishes

Torsten.

Answer 4

[Aldo]

Aren't you suppose to use this? Where did you get x_(i+1/2)=(x_(i+1)+x_(i))/2 x_(i-1/2)=(x_(i)+x_(i-1))/2 from?

Best regards Aldo :D

Answer 5

[Aldo]

Can't figure out a way to get the equation to look like this. And shouldn't (dT/dx)(x+h/2) ~ (T(x+h)-T(x))/h, be divided by 2h?

Best regards

Answer 6

[Aldo]

 clear all; close all, clf;

 T0=300; Tslut=410;
 N=100;
 h=410/N; n=N-1; x=h*(0:n)'; X=[0;x;410];
 k=2+x/7; Q=260*exp(-(x-410/2).^2); 

 A= sparse(n);
 for i=1:n-1, A(i,i)=(2*k/h^2); end
 for i=1:n-1, A(i,i+1)= (k/h^2); A(i+1,i)= (k/h^2);end 

 b= Q
 (1)= b(1)+(k/h^2)*T0
 b(n)= b(n)+(k/h^2)*Tslut

 T= A/b

I get this error "Subscripted assignment dimension mismatch. Error in uppgift5 (line 9) for i=1:n-1, A(i,i)=(2*k/h^2); end "

can you see how I can resolve it?

[Torsten]

Try

L=3.4;
N=100;
T0=300;
Tslut=410;
h=L/N;
fk=@(x)2+x/7;
fQ=@(x)260*exp(-(x-L/2)^2); 
A=zeros(N+1);
b=zeros(N+1,1);
A(1,1) = 1.0;
b(1) = T0;
for i=2:N
  x=h*(i-1);
  xhl=x-h/2;
  xhr=x+h/2;
  A(i,i-1) = -fk(xhl)/h^2; 
  A(i,i)   = (fk(xhl)+fk(xhr))/h^2;
  A(i,i+1) = -fk(xhr)/h^2;
  b(i) = fQ(x);
end
A(N+1,N+1) = 1.0;
b(N+1) = Tslut;
T=A\b;

Best wishes

Torsten.

[Aldo]

Hi,

isn't it better to use sparse(n)? Otherwise it takes ages to compute the result. I am trying to find the Temperature T at x = 1.5 with 6 decimals, and solve the equation. How would you plot the graph?

Best regards Aldo :D

[Torsten]

L=3.4;
N=100;
T0=300;
Tslut=410;
h=L/N;
x=linspace(0,L,N+1);
x=x.';
fk=@(x)2+x/7;
fQ=@(x)260*exp(-(x-L/2)^2); 
A=zeros(N+1);
b=zeros(N+1,1);
A(1,1) = 1.0;
b(1) = T0;
for i=2:N
  xp=x(i);
  xhl=xp-h/2;
  xhr=xp+h/2;
  A(i,i-1) = -fk(xhl); 
  A(i,i)   = (fk(xhl)+fk(xhr));
  A(i,i+1) = -fk(xhr);
  b(i) = fQ(xp)*h^2;
end
A(N+1,N+1) = 1.0;
b(N+1) = Tslut;
T=A\b;
plot(x,T)

If it really takes so much time to solve the linear system (I doubt it), I leave it to you to create a sparse version.

Best wishes

Torsten.

[Aldo]

Yes you are right it doesn't take much time to calculate. However, how do you get T specifically at x=1.5? I am interested in finding the value of T at point 1.5

Best Regards Aldo :D

[Torsten]

Take N=340 - then T at 1.5 is T(151) :-)

For arbitrary N,

T15 = interp1(x,T,1.5);

Best wishes

Torsten.