1-D Convection Equation is converging to 0. Thomas Algoritm using Back Time Center Scheme

1 次查看(过去 30 天)
Hi guys, this is a very mathermatical question i guess.
I have the following code for plotting 1-D convection equation in time given boundary conditions W1 and W2.
For some reason that i don't see the solution is converging to 0.
Please let me know if you see any errors if you are a MATH person
%Constants
L=10; %Lenght of interval x
a=0;
b=L;
alpha=0.2 %Diffusion
v=3; %Velocity
Nx=100; %Number of intervals
Tmax=10; %Max time
dt=0.01; %time step
Nt=Tmax/dt;
My_Fig=500;
plot_every=10;
%% Setup
dx=(b-a)/Nx;
xv=a:dx:b;
beta=alpha*dt/(dt^2);
gamma=v*dt/(2*dx);
u_ini=ff(xv);% Initial Condition
u_now=u_ini; %current time slice
plot_solution(xv,u_ini,u_now,My_Fig,0,0);
u_next=zeros(size(u_now));
A=zeros(Nx-1,Nx-1);
for i=1:Nx-1
A(i,i)=1+2*beta;
A(i,i+1)=-gamma-beta;
A(i+1,i)=gamma-beta;
end
A=A(1:Nx-1,1:Nx-1);
%% Iteration
Begin=now;
u_ini2=u_now;
for n=1:Nt
u_next(1)=u_now(1)-gamma*(u_now(2)-u_now(end-1))+beta*(u_now(2)-2*u_now(1)+u_now(end-1));
u_next(end)=u_now(1)-gamma*(u_now(2)-u_now(end-1))+beta*(u_now(2)-2*u_now(1)+u_now(end-1));
bvec=u_now(2:end-1)';
bvec(1)=bvec(1)-gamma*u_next(1)+beta*u_next(1);
bvec(end)=bvec(end)+gamma*u_next(end)+beta*u_next(end);
solvevec=inv(A)*bvec;
u_now=u_next;
if mod(n,plot_every)==0
plot_solution (xv,u_ini,u_now,My_Fig,0,dt*n);
end
end
End=now;
fprintf('BTCS: dt=%2.e Nx=%d/n',dt,Nx)
function w = ff(x)
% constructed so f(0)=f(10) and f'(0)=f'(10)=0
w1 = ((10/(3*pi))*sin(x*3*pi/10)+(0.1*(x-5).^2));
w2 = w1;
% w =(x-5).^2 / 10 - 10/3*sin(x*2*pi*1.5/10);
w = w1+w2-5;
return
end
function plot_solution(xv,u1,u2,FigNo,HOLD,time_now)
if isempty(FigNo), FigNo=1; end
if isempty(HOLD), HOLD = 0; end
figure(FigNo);
if HOLD, hold on; end
plot(xv,u1,'r-',xv,u2,'bo-');
if ~isempty(time_now)
tstring = sprintf('Thomas: T=%f',time_now);
title(tstring);
xlabel 'x';
ylabel 'u'(x,t);
end
%FormatThisFigure;
drawnow;
return
end

回答(1 个)

darova
darova 2020-5-4
Try this solution

类别

Help CenterFile Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by