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.
采纳的回答
Torsten
2016-11-8
(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.
更多回答（5 个）
Torsten
2016-11-8
[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.
Torsten
2016-11-9
-(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.
Aldo
2016-11-8
0 个投票
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
Aldo
2016-11-8
0 个投票
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
Aldo
2016-11-9
5 个评论
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
2016-11-9
Torsten
2016-11-10
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.
Torsten
2016-11-10
Take N=340 - then T at 1.5 is T(151) :-)
For arbitrary N,
T15 = interp1(x,T,1.5);
Best wishes
Torsten.
类别
在 帮助中心 和 File Exchange 中查找有关 Geometric Transformation and Image Registration 的更多信息
Translated by 
