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

1 回表示 (過去 30 日間)
Diego Dranuta
Diego Dranuta 2020 年 5 月 4 日
回答済み: darova 2020 年 5 月 4 日
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 Center および File ExchangeProgramming についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by