%global x_vec dx k
%clear
%clc
%clear all
%%%% du/dt = d^2u/dx^2; u is f(x,t)
%%%% DFM: (u(x, t+dt) -u(x, t))/dt = k*(u(x, t+dt) -2*u(x, t) - u(x-dt, t))/dx^2
%%%% Conctant k
k =2;
%%%% Length of Pipe
L= 10;
%%%% Number of elements
M = 10;
%%%% Discretize xspace
x_vec = linspace(0, L, M);
dx = x_vec(2) - x_vec(1);
%%%% Discretize t
dt = 0.5*(dx^2)/(2*k);
t_vec = 0:dt:10;
%%%% Allocate memory for u
u_mat = zeros(length(x_vec), length(t_vec));
%%%%% IC
u_mat(1, :) =200; % Left end of pipe
u_mat(end, :) = 150; % right end of the pipe
%%%% Euler and FM
for tdx = 1:length(t_vec)-1 %%%% integral inteval
for idx = 2:length(x_vec)-1
u_mat(idx, tdx+1) = u_mat(idx, tdx) + k*dt/(dx^2)*(u_mat(idx+1, tdx) -2*u_mat(idx, tdx) + u_mat(idx-1, tdx));
end
end
%%%% Plot
[tt, xx] = meshgrid(t_vec, x_vec);
mesh(xx, tt, u_mat)
xlabel('x')
ylabel('time')
zlabel('u')
%%%% Discretize t
dtrk4 = 0.5*(dx^2)/(2*k);
trk4_vec = t_vec(1): dtrk4:t_vec(end);
% Allocate memory for u
urk4_mat = zeros(length(x_vec), length(trk4_vec));
%%%% IC
urk4_mat(1,:) = 200; % Left end of pipe
urk4_mat(end,:) = 150; % right end of the pipe
for tdx = 1:length(trk4_vec)-1
k1 = Derivative(urk4_mat(:, tdx));
k2 = Derivative(urk4_mat(:, tdx) + k1*dtrk4/2);
k3 = Derivative(urk4_mat(:, tdx) + k2*dtrk4/2);
k4 = Derivative(urk4_mat(:, tdx) + k3*dtrk4);
phi = (1/6)*(k1+ 2*k2 + 2*k3 +k4);
urk4_mat(:, tdx+1) = urk4_mat(:, tdx) + phi*dtrk4;
end
%%%% Plot this
[tt, xx] = meshgrid(trk4_vec, x_vec);
figure()
mesh(xx, tt, urk4_mat)
xlabel('x')
ylabel('time')
zlabel('u')
function dudt = Derivative(Tin_vec)
%global x_vec dx k
dudt = 0*Tin_vec;
x_vec = linspace(0,10,10);
k = 2;
dx = x_vec(2) - x_vec(1);
for idx = 2:length(x_vec)-1
dudt(idx) = k/(dx^2)*(Tin_vec(idx+1) -2*Tin_vec(idx) + Tin_vec(idx-1));
end
end
