burgers_exact.m
代码语言：javascript复制function burgers_exact
% Plot the exact solution of Burgers' equation
% u_t u*u_x = 0, u(x,0)=f(x)
u = 0.01:0.01:0.999;
xp = @(u,t) u.*t sqrt((1 - u)./u) ;
xm = @(u,t) u.*t - sqrt((1 - u)./u) ;
plot(xp(u,0),u,'-b',xm(u,0),u,'-b','LineWidth',2), grid on
xlabel 'x', ylabel 'u'
title ' Solution of u_{t} uu_{x} = 0'
for t = 0:1:10
plot(xp(u,t),u,'-k',xm(u,t),u,'-k','LineWidth',2), grid on
pause(0.2);
end
burgers_exact2.m
代码语言：javascript复制function burgers_exact2
%SIDENOTE
%If the characteristics of the exact solution do not intersect, then the
%classical solution of burgers equation is given by:
% Domain
x = 0:0.02:1;
% IC (at t=0)
u0 = 0.2 sin(2*pi*x);
% u function
u = @(x,t) 0.2 sin(2*pi*x);
for t = 0:0.01:0.1
u_next = u(x - t*u(x,t),t);
plot(x,u_next,'o');
pause(0.1)
drawnow
end
% This solution is valid for any u(x,t) function.
% However if the characteristics colide, then the classical solution is not
% valid anymore. It can only be solved computationally.LinearAdvection_1D.m
代码语言：javascript复制%% Linear 1D Advection equation
% Routine for solving the following PDE
%
% du/dt a du/dx = 0
%
% Using the following scheemes:
%
% 1. One-sided Upwind
% 2. Lax-Friedrichs
% 3. Lax-Wendroff
% 4. Beam-Warming
%
clc; clear; close all;
%% Parameters
cfl = 0.8; % CFL = a*dt/dx
tend = 0.2; % End time
a = 0.5; % Scalar wave speed
%% Domain Discretization
dx = 0.01;
dt = cfl/abs(a)*dx;
x = 0:dx:1;
t = 0:dt:tend;
%% Initial Condition
n = length(x);
u_0 = zeros(1,n);
u_0(1:floor(n/2)) = 1; u_0(floor(n/2) 1:n) = 0;
u_next = zeros(1,n);
%% 1. One-Sided Upwind
u_upwind = u_0;
for k = t
for j = 2:n-1
u_next(j) = u_upwind(j) - cfl*(u_upwind(j)-u_upwind(j-1));
end
u_upwind = u_next;
% BC
u_upwind(1) = u_next(2);
u_upwind(n) = u_next(n-1);
end
u_next = zeros(1,n);
%% 2. Lax-Friedrichs
u_LF = u_0;
for k = t
for j = 2:n-1
u_next(j) = 1/2*(u_LF(j-1) u_LF(j 1)) - 1/2*cfl*(u_LF(j 1)-u_LF(j-1));
end
u_LF = u_next;
% BC
u_LF(1) = u_next(2);
u_LF(n) = u_next(n-1);
end
u_next = zeros(1,n);
%% 3. Lax-Wendroff
u_LW = u_0;
for k = t
for j = 2:n-1
u_next(j) = u_LW(j) - 1/2*cfl*(u_LW(j 1)-u_LW(j-1)) ...
1/2*(cfl^2)*(u_LW(j 1)-2*u_LW(j) u_LW(j-1));
end
u_LW = u_next;
% BC
u_LW(1) = u_next(2);
u_LW(n) = u_next(n-1);
end
u_next = zeros(1,n);
%% 4. Beam-Warming
u_BW = u_0;
for k = t
for j = 3:n
u_next(j) = u_BW(j) - 1/2*cfl*(3*u_BW(j)-4*u_BW(j-1) u_BW(j-2)) ...
1/2*(cfl^2)*(u_BW(j)-2*u_BW(j-1) u_BW(j-2));
end
u_BW = u_next;
% BC
u_BW(1) = u_next(3);
u_BW(2) = u_next(3);
end
%% Exact Solution
x_exact = 0.5 a*tend;
u_exact = zeros(1,n);
for j = 1:n
if x(j) <= x_exact
u_exact(j) = 1;
else
u_exact(j) = 0;
end
end
%% Make a comparative Plot
offset = 0.05;
subplot(221);
hold on
plot(x,u_upwind,'.'); plot(x,u_exact,'k'); xlabel('X-Coordinate [-]'); ylabel('U-state [-]'); ylim([-0.5,1.5]); title 'One-Sided Upwind';
hold off
subplot(222);
hold on
plot(x,u_LF,'.'); plot(x,u_exact,'k'); xlabel('X-Coordinate [-]'); ylabel('U-state [-]'); ylim([-0.5,1.5]); title 'Lax-Friedrichs';
hold off
subplot(223);
hold on
plot(x,u_LW,'.'); plot(x,u_exact,'k'); xlabel('X-Coordinate [-]'); ylabel('U-state [-]'); ylim([-0.5,1.5]); title 'Lax-wendroff';
hold off
subplot(224);
hold on
plot(x,u_BW,'.'); plot(x,u_exact,'k'); xlabel('X-Coordinate [-]'); ylabel('U-state [-]'); ylim([-0.5,1.5]); title 'Beam-Warming';
hold offsweby.m
代码语言：javascript复制%% Plot sweby diagrams
% number of points
N = 1000;
% compute flux limiting functions for a range of methods
theta = linspace(-1, 5, N);
upwind = zeros(N,1);
laxwendroff = ones(N,1);
beamwarming = theta;
fromm = 0.5*(1 theta);
minmo = minmod(1, theta);
superbee = max(0, max(min(1,2*theta), min(2,theta)));
MC = max(0, min(min( (1 theta)/2, 2), 2*theta));
vanleer = (theta abs(theta))./(1 abs(theta));
% plot TVD stability region
% vertices of polygonal bounding region (truncated at x=10)
figure(1)
x = [0 10 10 1];
y = [0 0 2 2];
fill(x,y,[.9 .9 .9]); % fill region with grey
hold on;
plot(theta, upwind, 'r-', 'LineWidth', 1.5)
plot(theta, laxwendroff, 'b-', 'LineWidth', 1.5)
plot(theta, beamwarming, 'g-', 'LineWidth', 1.5)
plot(theta, fromm, 'm', 'LineWidth', 1.5)
hold off;
axis equal
axis([-1 4 -.5 2.5])
legend('TVD region', 'Upwind', 'Lax-Wendroff', 'Beam-Warming', 'Fromm');
xlabel('theta', 'FontSize', 18);
ylabel('phi', 'FontSize', 18);
% plot Sweby modified TVD stability region
% vertices of polygonal bounding region (truncated at x=10)
figure(2)
x = [0 1 10 10 2 1 0.5 0];
y = [0 1 1 2 2 1 1 0];
fill(x,y,[.9 .9 .9]); % fill Sweby TVD region with grey
hold on;
plot(theta, minmo, 'r-', 'LineWidth', 1.5);
plot(theta, superbee, 'g-', 'LineWidth', 1.5);
plot(theta, MC, 'b-', 'LineWidth', 1.5);
plot(theta, vanleer, 'm-', 'LineWidth', 1.5);
hold off;
axis equal
axis([-1 4 -.5 2.5])
legend('Sweby TVD region', 'minmod', 'superbee', 'MC', 'van Leer');
xlabel('theta', 'FontSize', 18);
ylabel('phi', 'FontSize', 18);
TVD_scheme.m
代码语言：javascript复制%% 1D Linear Advection Equation
%
% du/dt a du/dx = 0
%
% Subroutine for solving using TVD Method.
% by Manuel Diaz, manuel.ade'at'gmail.com
% Institute of Applied Mechanics, 2012.08.12
clear all; close all; clc;
%% Parameters
a = -0.5;
a_m = min(0,a);
a_p = max(0,a);
dx = 0.01;
cfl = 0.9;
dt = cfl*dx/abs(a);
t_end = 0.4;
%% Discretization of Domain
x = 1:dx:2;
t = 0:dt:t_end;
%% Initial Condition
n = length(x);
%u_0 = zeros(1,n);
u_0(1:ceil(n/2)) = 1;
u_0((ceil(n/2) 1):n) = 2;
%% Main Loop
% Initilize vector variables
r = zeros(1,n);
F_rl = zeros(1,n);
F_rh = zeros(1,n);
F_ll = zeros(1,n);
F_lh = zeros(1,n);
F_right = zeros(1,n);
F_left = zeros(1,n);
u_next = zeros(1,n);
% Load Initial condition
u = u_0;
% Start TVD Loop
for k = t
for j = 2:n-1
% smooth measurement factor 'r'
if u(j) == u(j 1)
r(j) = 1;
elseif a > 0
r(j) = (u(j) - u(j-1)) / (u(j 1) - u(j));
elseif a < 0
r(j) = (u(j 2) - u(j 1)) / (u(j 1) - u(j));
end
r(1) = 1; r(n) = 1;
end
% Define Flux Limiter function
phi = (r abs(r))./(1 abs(r));
for j = 2:n-1
% Compute fluxes for TVD
F_rl(j) = a_p*u(j) a_m*u(j 1);
F_rh(j) = (1/2)*a*(u(j) u(j 1)) - (1/2)*(a^2)*(dt/dx)*(u(j 1)-u(j));
F_ll(j) = a_p*u(j-1) a_m*u(j);
F_lh(j) = (1/2)*a*(u(j-1) u(j)) - (1/2)*(a^2)*(dt/dx)*(u(j)-u(j-1));
% Compute next time step
F_right(j) = F_rl(j) phi(j)*( F_rh(j) - F_rl(j) );
F_left(j) = F_ll(j) phi(j-1)*( F_lh(j) - F_ll(j) );
u_next(j) = u(j) - dt/dx*(F_right(j) - F_left(j));
end
% right sided Upwind
%u_next(j) = u(j) - a*dt/dx*(u(j)-u(j-1));
% letf sided Upwind
%u_next(j) = u(j) - a*dt/dx*(u(j 1)-u(j));
% Lax-Wendroff
%u_next(j) = u(j) - dt/(2*dx)*a*(u(j 1)-u(j-1)) ...
% (dt^2/2)/(dx^2)*(a^2)*(u(j 1)-2*u(j) u(j-1)); % LW
% Double sided Upwind
%u_next(j) = u(j)-a_p*dt/dx*(u(j)-u(j-1))-a_m*dt/dx*(u(j 1)-u(j));
% BC
u_next(1) = u_next(2);
u_next(n) = u_next(n-1);
% UPDATE info
u = u_next(1:n);
end
%% Plot Results
plot(x,u,'.'); Ylim([0.5,2.5]); Title('TVD');
xlabel('X-Coordinate [-]'); ylabel('U-state [-]');Upwind.m
代码语言：javascript复制%% 1D Linear Advection Equation
%
% du/dt a du/dx = 0
%
% Ref's:
% Randall J. Leveque;
% Finite Volume Methods for Hyperbolic Problems
% Cambridge Press, 2002.
%
% Subroutine for solving using Upwind Method.
% by Manuel Diaz, manuel.ade'at'gmail.com
% Institute of Applied Mechanics, 2012.08.12
clear all; close all; clc;
%% Parameters
a = 0.5;
a_m = min(0,a);
a_p = max(0,a);
dx = 0.01;
cfl = 0.9;
dt = cfl*dx/abs(a);
t_end = 0.6;
%% Discretization of Domain
x = 1:dx:2;
t = 0:dt:t_end;
%% Initial Condition
n = length(x);
%u_0 = zeros(1,n);
u0(1:ceil(n/2)) = 1;
u0((ceil(n/2) 1):n) = 2;
%% Main Loop
% Initilize vector variables
u_next = zeros(size(u0)); % u^{n 1}
u_bar1 = zeros(size(u0)); % preallocate for MacCormack
u_bar2 = zeros(size(u0)); % preallocate for MacCormack
% Load Initial condition
u = u0;
% Main Loop
for k = t
% plot every time step
plot(x,u,'.')
for j = 2:n-1
% Single Sided Upwind
%u_next(j) = u(j) - a*dt/dx*(u(j)-u(j-1)); % Upwind
%u_next(j) = u(j) - a*dt/dx*(u(j 1)-u(j)); % Downwind
% Double Sided Upwind
u_next(j) = u(j) - a_p*dt/dx*(u(j)-u(j-1)) - a_m*dt/dx*(u(j 1)-u(j));
% Comparing to Lax-Wendroff:
%u_next(j) = u(j) - a*dt/(2*dx)*(u(j 1)-u(j-1)) ...
% (dt^2/2)/(dx^2)*(a^2)*(u(j 1)-2*u(j) u(j-1)); % LW
% Comparing to MacCormack Method:
% Predictor steps
%u_bar1(j) = u(j) - a*dt/dx*(u(j 1)-u(j));
% u_bar1(1) = u_bar1(2); % Neumann BC in predictor step
%u_bar2(j) = u_bar1(j) - a*dt/dx*(u_bar1(j)-u_bar1(j-1));
% Corrector step
%u_next(j) = 1/2*(u(j) u_bar2(j));
end
% BC
u_next(1) = u_next(2); % Neumann BC
u_next(n) = u_next(n-1); % Neumann BC
% UPDATE info
u = u_next(1:n);
% update figure
pause(0.05); drawnow
end
%% Plot final Result
plot(x,u,'.')Upwind_for_Nazanin.m
代码语言：javascript复制%% 1D Linear Advection Equation
% for Nazanin
%
% du/dt a du/dx = 0
%
% Subroutine for solving using Upwind Method.
% by Manuel Diaz, manuel.ade'at'gmail.com
% Institute of Applied Mechanics, 2013.03.08
clear all; close all; clc;
%% Parameters
dx = 1/200; % using 200 points
cfl = 0.9;
t_end = 0.5;
%% Discretization of Domain
x = 0:dx:1;
n = length(x);
%% Initial Condition
u0 = zeros(1,n);
idx = find( x>=0.2 & x<=0.4);
u0(idx) = 1;
%% Main Loop
% Initilize vector variables
u_next = zeros(size(u0)); % u^{n 1}
% Load Initial condition
u = u0;
% Main Loop
t = 0;
while t < t_end
% plot every time step
plot(x,u,'.')
% Compute advection speed
a = (1 x.^2)./(1 2*x*t 2*x.^2 x.^4);
a_m = min(0,a);
a_p = max(0,a);
% compute time step
dt = cfl*dx/abs(max(a));
t = t dt;
for j = 2:n-1
% Using double Sided Upwind
u_next(j) = u(j) - a_p(j)*dt/dx*(u(j)-u(j-1)) - a_m(j)*dt/dx*(u(j 1)-u(j));
end
% BC
u_next(1) = 0; % Dirichlet BC
%u_next(n) = u_next(n-1); % Neumann BC
% UPDATE info
u = u_next;
% update figure
pause(0.01); drawnow
end
%% Compute exact solution at t = t_end
u_exact = zeros(1,n);
idx = find( x>=(0.2 - t_end/(1 0.2^2) ) & x<=(0.4 - t_end/(1 0.2^2)));
u_exact(idx) = 1;
%% Plot final Result
hold on
plot(x,u,'.')
%plot(x,u_exact,'-k')
hold off
http://mpvideo.qpic.cn/0b78leaacaaaymac5mocazpfawodafmqaaia.f10002.mp4?dis_k=1828ad8a9e8303d0d69ab2c8ca6bfd22&dis_t=1653729006&vid=wxv_1345401917991862273&format_id=10002&support_redirect=0&mmversion=false
