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