I need help solving a PDE using Thomas Method

8 views (last 30 days)
Hi,
I am trying to use MATLAB for solving PDEs for a MATH course I am taking but I am having issues with the coding section because I am not very familiar with MATLAB.
I attached the problem and the code I have developed so far. The issue I am havinf is to plot the solution that varies with time. I have solve the math part but plotting is not my strenght. Can you take a look and tell me how I can proceed?
Thank you!
and here is my code so far
%Constants
L=10; %Lenght of interval x
a=0;
b=L;
alpha=0.2 %Diffusion
v=3; %Velocity
Nx=200; %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=fft(xv);% Initial Condition
u_now=u_ini; %current time slice
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_uni2=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;
end
End=now;
fprintf('BTCS: dt=%2.e Nx=%d/n',dt,Nx)

Answers (0)

Tags

Community Treasure Hunt

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

Start Hunting!