function linearAdvection_JAclear % Clears workspace
clc % Clears history
clf % Clears figure window
N = 101; % Defines the no. of grid points
p = 0; % Sets the lower end of the domain
q = 100; % Sets the upper end of the domain
v = 0.8; % defines velocity of the water (meters)
dx = (q - p)/(N - 1); % Calculates the grid spacing
dt = 0.8; % Defines the time step.
t = 0; % Sets the initial time to zero
runtime = 100; % Defines time program runs for
s = 0.4; % Sets safety factor to 0.4 (Advised not to change this)
pausetime = 0.1; % Pause time between animation
L=1;x = p : dx : q; % Calculates the value of x using the domain and the grid spacing
for i = 1 : N U(i) = exp(-(x(i)-L/2).^2); % Defines the initial pollutant using initial conditions
end% fprintf('Select a scheme from the list and enter the corrosponding number \n\nFirst Order Upwind - 1\nLax-Friedrichs - 2 \nLax-Wendroff - 3\nAdams Average (Lax-Friedrichs) - 4\nAdams Average (Lax-Wendroff) - 5\n\n')
Choice = input('Scheme no. = '); % Select a scheme to run
if Choice == 2 A = 1; B = 1; % <<<< DO NOT CHANGE THESE VALUES
C = 0;elseif Choice == 3 A = 1; B = 1; % <<<< DO NOT CHANGE THESE VALUES C = 0;elseif Choice == 4 A = 2; B = 2; % <<<<< Values can be changed
C = 1;elseif Choice == 5 A = 1; B = 1; % <<<<<<<< Values can be changed
C = 5;end% ===================== Runs similation ===================================
while t < runtime % Halts program when inequality is violated
% Applies periodic boundary conditions depending on choice
if Choice == 1 boundary1 = un(N); boundary2 = un(1); un = [boundary1 un(1:N-1) boundary2]; elseif Choice == 2 boundary1 = un(N); boundary2 = un(1); un = [boundary1 un boundary2]; elseif Choice == 3 boundary1 = un(N); boundary2 = un(1); un = [boundary1 un boundary2]; elseif Choice == 4 boundary1 = un(N); boundary2 = un(1); un = [boundary1 un boundary2]; elseif Choice == 5 boundary1 = un(N); boundary2 = un(1); un = [boundary1 un boundary2]; end dt = s * dx / v; % Calculates time step
t = t + dt; % Adds time step to time
c = dt/dx; % Calculates this here to provide efficiency
if Choice == 1 % Applies First Order Upwind scheme
for j = 2 : N + 1 u(j) = un(j) - v * c * (un(j) - un(j-1)); end elseif Choice == 2 % Applies either Lax-Friedrichs or Adams Average (Lax-Friedrichs) scheme
for j = 2 : N + 1 u(j) = (A*un(j+1) + B*un(j-1) + C*un(j))/(A+B+C) - v * c * (un(j+1) - un(j-1)); end elseif Choice == 4 % Applies either Lax-Friedrichs or Adams Average (Lax-Friedrichs) scheme for j = 2 : N + 1 u(j) = (A*un(j+1) + B*un(j-1) + C*un(j))/(A+B+C) - v * c * (un(j+1) - un(j-1)); end elseif Choice == 3 % Applies either Lax-Wendroff or Adams Average (Lax-Wendroff) scheme
for j = 2 : N + 1 u(j) = (A*un(j+1) + B*un(j-1) + C*un(j))/(A+B+C) - v * c * (un(j+1) - un(j-1))/2 + c^2 * v^2 * 0.5 * (un(j+1) - 2*un(j) + un(j-1)); end elseif Choice == 5 % Applies either Lax-Wendroff or Adams Average (Lax-Wendroff) scheme for j = 2 : N + 1 u(j) = (A*un(j+1) + B*un(j-1) + C*un(j))/(A+B+C) - v * c * (un(j+1) - un(j-1))/2 + c^2 * v^2 * 0.5 * (un(j+1) - 2*un(j) + un(j-1)); end end plot(x,u(2:N+1),'b') % Plots pollutant (dependent variable)
xlabel('x [m]') % Adds appropiate labels
ylabel('U [m]') % Adds title depending on scheme choice
if Choice == 1 title('Solves 1D Advection Equation using FOU Scheme','Fontsize',10) elseif Choice == 2 title('Solves 1D Advection Equation using Lax-Friedrichs Scheme','Fontsize',10) elseif Choice == 3 title('Solves 1D Advection Equation using Lax-Wendroff Scheme','Fontsize',10) elseif Choice == 4 title('Solves 1D Advection Equation using Adams Average (Lax-Friedrichs) Scheme','Fontsize',10) elseif Choice == 5 title('Solves 1D Advection Equation using Adams Average (Lax-Wendroff) Scheme','Fontsize',10) end axis([p q min(u) 1.1*max(u)]) % Sets axis
pause(pausetime) % Shows results for each timestep
un = u(2:N+1); % Ommits boundary conditions ready for the next time step
endendfunction initial = U(x) % Defines initial conditions.
pollutantstart1 = 20; % Start of pollutant concentration 1
pollutantend1 = 40; % End of pollutant concentration 1
pollutantstart2 = 60; % Start of pollutant concentration 2
pollutantend2 = 80; % End of pollutant concentration 2
pollutantfloor = 0.1; % Defines bottom of pollutant
pollutantpeak = 0.8; % Defines top of pollutant
initial = pollutantfloor; % Sets up bottom of pollutant
% ============ Determins the areas in which the pollutant peaks ==========
if x >= pollutantstart1 && x <= pollutantend1 initial = pollutantpeak;elseif x >= pollutantstart2 && x <= pollutantend2 initial = pollutantpeak;endend
MATLAB: Error when i want to run it
please help me work it
Best Answer