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

2 views (last 30 days)
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

Categories

Find more on Get Started with MATLAB in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!