I tried Euler's method and RK4 method but I'm getting different results. Please can someone help
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponses (1)
Torsten
le 3 Août 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
Voir également
Catégories
En savoir plus sur Geometry and Mesh dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!