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

1 view (last 30 days)
Diego Dranuta
Diego Dranuta on 4 May 2020
Answered: darova on 4 May 2020
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

Answers (1)

darova
darova on 4 May 2020
Try this solution

Tags

Community Treasure Hunt

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

Start Hunting!