I'm writting a code to solve the "equation of advection", which express how a given property or physical quantity varies with time. For that I've this expression: https://gyazo.com/f531cb756ffbd3ec28ab85ea1f09b18d This is one of my issues, I dont know how to implement it. This is the code I've for the moment:
%%Our paramatres
k = 1e-4;
x = [1 2 2.1 13.9 14 28 28.1 39.9 40 59 69 79 150];
y = [3000 3000 3150 3150 3000 3000 3080 3080 3000 3000 3150 3000 3000];
vx = 0.1; % velocity in x m/s
tmax = 250; %time max
xmin = min(x);
xmax = max(x);
axis([0 150 2990 3160])
plot(x,y,'--')
%% discretization of domain
%
dx = 150;
dt = dx/(5*vx); % delta time is equal delta x * vx
n = tmax / dt ;
m = (xmax - xmin)/dx;
x = linspace(xmin,xmax,m+1); % é igual a x = xmin:dx:xmax
%% Stability condition
%
lambda = k * dt/vx^2;
while lambda > 1/2
fprintf ('Satisfied\n\n')
fprintf ('Value of dt is %5.1f and for z is %i \n\n', dt,dx)
dt = input('New value for dt: ')
dx = input ('New value for dx: ')
lambda = k*dt/dx^2;
n = tmax/dt;
m = xmax/dx;
x = linspace(xmin,xmax,m+1);
end
%%Initial conditions and frontier
%
h = zeros(size(x));
h(x>0)=1;
plot(x,h)
xlabel('Distance')
ylabel('Density')
grid on
