first of all I'd like to point out that this is my first time using Stack Overflow so don't hesitate to tell me if something isn't displaying properly or isn't clear.
I'm trying to solve an advection equation with a second member whose analytical solution is known. I also wanted to test the method with a finite difference discretization. The results are similar, but there is a difference. I calculate the solution over four portions. For the first portion, the exact solution coincides with the analytical solution. For the second portion, the agreement is quite good, but errors still appear. The situation deteriorates for the third and fourth portions.
I've already asked this question on stackoverflow and they advised me to look at this platform because from a code point of view there's no problem.
Here is the function that calculates the solutions analytically
function [C1, C2, C3, C4] = calcul_sol_analy(mesures,params)
% Solution portion 1
x_p1 = linspace(0, params.xcm, 100);
[X1, T1] = meshgrid(x_p1, params.t_values);
C1 = params.c0(T1 - X1/params.v) .* exp(-params.k * X1 / params.v);
% Solution portion 2
Cm = interp1(params.date_jours_ch, mesures, params.t_values, 'pchip')'; % transposée
x_p2 = linspace(params.xcm, params.xcp, 100);
[X2, T2] = meshgrid(x_p2, params.t_values);
c1 = @(t) (C1(:, end) .* params.debit_s + Cm .* params.debit_m) ./ (params.debit_SM);
C2 = c1(T2 - (X2-params.xcm)/params.v) .* exp(-params.k * (X2-params.xcm)/ params.v);
% Solution portion 3
Cp = interp1(params.date_jours_cl, params.mesures_crassphage_ch, params.t_values, 'pchip')'; % transposée
x_p3 = linspace(params.xcp, params.xco, 100);
[X3, T3] = meshgrid(x_p3, params.t_values);
c2 = @(t) (C2(:, end) .* params.debit_SM + Cp .* params.debit_cl) ./ (params.debit_SME);
C3 = c2(T3 - (X3-params.xcp)/params.v) .* exp(-params.k * (X3-params.xcp)/ params.v);
% Solution portion 4
Co = interp1(params.date_jours_tr, params.mesures_crassphage_ch, params.t_values, 'pchip')'; % transposée
x_p4 = linspace(params.xco, params.xct, 100);
[X4, T4] = meshgrid(x_p4, params.t_values);
c3 = @(t) (C3(:, end) .* params.debit_SME + Co .* params.debit_o) ./ (params.debit_SMO);
C4 = c3(T4 - (X4-params.xco)/params.v) .* exp(-params.k * (X4-params.xco)/ params.v);
end
and here is the function that calculates the solutions numerically
function [C1, C2, C3, C4] = calcul_sol_num(mesures, params)
% Solution portion 1
C1 = zeros(params.Nt, params.Nx);
C1(:, 1) = params.c0(params.t_values);
x_p1 = linspace(0, params.xcm, params.Nx);
[X1, T1] = meshgrid(x_p1, params.t_values);
C_exact = params.c0(T1 - X1 / params.v) .* exp(-params.k * X1 / params.v);
C1(1, :) = C_exact(1, :);
for n = 1:params.Nt-1
for i = 2:params.Nx-1
C1(n + 1, i) = C1(n, i) * (1 - params.k * params.dt - params.v * params.dt / params.dx) + ...
C1(n, i-1) * (params.v * params.dt / params.dx);
end
C1(n + 1, params.Nx) = C1(n + 1, params.Nx-1);
end
% Solution portion 2
Cm = interp1(params.date_jours_ch, mesures, params.t_values, 'pchip')'; % transposée
x_p2 = linspace(params.xcm, params.xcp, params.Nx);
[X2, T2] = meshgrid(x_p2, params.t_values);
c1 = @(t) (C1(:, end) .* params.debit_s + Cm .* params.debit_m) ./ params.debit_SM;
C2 = zeros(params.Nt, params.Nx);
C2(:, 1) = c1(params.t_values);
C_exact = c1(T2 - (X2-params.xcm) / params.v) .* exp(-params.k * (X2-params.xcm) / params.v);
C2(1, :) = C_exact(1, :);
for n = 1:params.Nt-1
for i = 2:params.Nx-1
C2(n + 1, i) = C2(n, i) * (1 - params.k * params.dt - params.v * params.dt / params.dx) + ...
C2(n, i-1) * (params.v * params.dt / params.dx);
end
C2(n + 1, params.Nx) = C2(n + 1, params.Nx-1);
end
% Solution portion 3
Cp = interp1(params.date_jours_cl, params.mesures_crassphage_ch, params.t_values, 'pchip')'; % transposée
x_p3 = linspace(params.xcp, params.xco, params.Nx);
[X3, T3] = meshgrid(x_p3, params.t_values);
c2 = @(t) (C2(:, end) .* params.debit_SM + Cp .* params.debit_cl) ./ params.debit_SME;
C3 = zeros(params.Nt, params.Nx);
C3(:, 1) = c2(params.t_values);
C_exact = c2(T3 - (X3-params.xcp) / params.v) .* exp(-params.k * (X3-params.xcp) / params.v);
C3(1, :) = C_exact(1, :);
for n = 1:params.Nt-1
for i = 2:params.Nx-1
C3(n + 1, i) = C3(n, i) * (1 - params.k * params.dt - params.v * params.dt / params.dx) + ...
C3(n, i-1) * (params.v * params.dt / params.dx);
end
C3(n + 1, params.Nx) = C3(n + 1, params.Nx-1);
end
% Solution portion 4
Co = interp1(params.date_jours_tr, params.mesures_crassphage_oi, params.t_values, 'pchip')';
x_p4 = linspace(params.xco, params.xct, params.Nx);
[X4, T4] = meshgrid(x_p4, params.t_values);
c3 = @(t) (C3(:, end) .* params.debit_SME + Co .* params.debit_o) ./ params.debit_SMO;
C4 = zeros(params.Nt, params.Nx);
C4(:, 1) = c3(params.t_values);
C_exact = c3(T4 - (X4-params.xco) / params.v) .* exp(-params.k * (X4-params.xco) / params.v);
C4(1, :) = C_exact(1, :);
for n = 1:params.Nt-1
for i = 2:params.Nx-1
C4(n + 1, i) = C4(n, i) * (1 - params.k * params.dt - params.v * params.dt / params.dx) + ...
C4(n, i-1) * (params.v * params.dt / params.dx);
end
C4(n + 1, params.Nx) = C4(n + 1, params.Nx-1);
end
end
Discretization parameter :
Nx = 101;
Tf = max(date_jours_ch);
dx = L / (Nx-1) ;
dt = dx / v;
x_values = 0:dx:L;
t_values = 2:dt:Tf;
Nt=length(t_values);
solution portion 1 :
solution portion 2 :
solution portion 3 :
solution portion 4 :
