线性一维平流方程解决以下PDE的程序。

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

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.

%% 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 off

%% 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);

%% 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 [-]');

%% 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,'.')

%% 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