I tried Euler's method and RK4 method but I'm getting different results. Please can someone help
1 view (last 30 days)
Show older comments
function HeatPDE()
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, 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;
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
0 Comments
Answers (1)
Torsten
on 3 Aug 2022
%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
See Also
Categories
Find more on Geometry and Mesh 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!