I need help solving a PDE using Thomas Method
8 views (last 30 days)
Show older comments
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)
0 Comments
Answers (0)
See Also
Categories
Find more on Boundary Conditions in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!