-

MATH_5343_Numerical_Solutions_to_PDE_Zeng/adv/adv_cauchy.m

@@ -0,0 +1,68 @@
+function [x, u] = adv_cauchy(N,cfl,c,dt,dx,fic)
+  % [x, u] = adv_cauchy(N,cfl,c,dt,dx,ic)
+  % Compute the solution to the Cauchy problem:
+  %    u_t + cu u_x = 0               0 <= x, t <= 1
+  %    u(0,t) = u(1,t)                0 <= t <= 1
+  %    u(x,0) is specified by the input "fic".
+  % using a uniform grid with N sub-intervals.
+  % Other arguments:
+  %  dt:  'be' | 'fe' | 'cn'
+  %  dxx: 'cd' | 'upw'
+
+  x = linspace(0,1,N+1);
+  h = 1./N; % spatial interval size
+  T = 1;
+  % Determine the time step size by the CFL condition 
+  k = cfl*h/abs(c);
+
+  % fic should be a function handle
+  u = fic(x');
+
+  % Spatial discretization
+  A = zeros(N+1);
+  I = eye(N+1);
+  % Setup index for periodic boundary condition
+  rght = @(i)(i<N+1)*(i+1)+(i==N+1)*2;
+  left = @(i)(i>1)*(i-1)+(i==1)*N;
+  % - discretization of advection part -
+  if strcmp( dx, 'cd' ) == 1
+    for i = 1 : N+1
+      A(i,left(i)) = A(i,left(i)) - .5*c/h;
+      A(i,rght(i)) = A(i,rght(i)) + .5*c/h;
+    end
+  elseif strcmp( dx, 'upw' ) == 1
+    if c<0
+      for i = 1 : N+1
+        A(i,i)       = A(i,i) - c/h;
+        A(i,rght(i)) = A(i,rght(i)) + c/h;
+      end
+    else
+      for i = 1 : N+1
+        A(i,left(i)) = A(i,left(i)) - c/h;
+        A(i,i)       = A(i,i) + c/h;
+      end
+    end
+  else
+    disp('The advection part is not supported');
+  end
+  % - function handle to march in time -
+  if strcmp( dt, 'fe' ) == 1
+    mt = @(kval,uval)(I-kval*A)*uval;
+  elseif strcmp( dt, 'be' ) == 1
+    mt = @(kval,uval)(I+kval*A)\uval;
+  elseif strcmp( dt, 'cn' ) == 1
+    mt = @(kval,uval)(I+.5*kval*A)\((I-.5*kval*A)*uval);
+  else
+    disp('The time-marching method is not supported');
+  end
+  % - solving the advection equation -
+  finish = 0;
+  while finish==0
+    if k>T % Check if the remaining time is smaller than k
+      k = T;
+      finish = 1;
+    end
+    u = mt(k,u); % Marching in time
+    T = T - k; % Remove time step size from remaining time.
+  end
+end

MATH_5343_Numerical_Solutions_to_PDE_Zeng/adv/adv_cauchy_lf.m

@@ -0,0 +1,43 @@
+function [x, u] = adv_cauchy_lf(N,cfl,c,fic)
+  % [x, u] = adv_cauchy_lf(N,cfl,c,ic)
+  % Compute the solution to the Cauchy problem:
+  %    u_t + cu u_x = 0               0 <= x, t <= 1
+  %    u(0,t) = u(1,t)                0 <= t <= 1
+  %    u(x,0) is specified by the input "fic".
+  % using a uniform grid with N sub-intervals and
+  % the Lax-Friedrich method.
+
+  x = linspace(0,1,N+1);
+  h = 1./N; % spatial interval size
+  T = 1;
+  % Determine the time step size by the CFL condition 
+  k = cfl*h/abs(c);
+
+  % fic should be a function handle
+  u = fic(x');
+
+  % Spatial discretization u^{n+1} = (D-kA) u^n
+  D = zeros(N+1);
+  A = zeros(N+1);
+  % Setup index for periodic boundary condition
+  rght = @(i)(i<N+1)*(i+1)+(i==N+1)*2;
+  left = @(i)(i>1)*(i-1)+(i==1)*N;
+  for i = 1 : N+1
+    D(i,left(i)) = .5;
+    D(i,rght(i)) = .5;
+  end
+  for i = 1 : N+1
+    A(i,left(i)) = A(i,left(i)) - .5*c/h;
+    A(i,rght(i)) = A(i,rght(i)) + .5*c/h;
+  end
+  % - solving the advection equation -
+  finish = 0;
+  while finish==0
+    if k>T % Check if the remaining time is smaller than k
+      k = T;
+      finish = 1;
+    end
+    u = (D-k*A)*u;
+    T = T - k; % Remove time step size from remaining time.
+  end
+end

MATH_5343_Numerical_Solutions_to_PDE_Zeng/adv/plot_example_cd_be.m

@@ -0,0 +1,32 @@
+c = 1.0;
+N = 80;
+cfl = {[0.99 2 5], [0.99 2 5]};
+cfl_spec = {'kx--','k*-.','ks:'};
+
+f_ic{1} = @(x)16*(x.*(1-x)).^2;
+f_ic{2} = @(x)(x<=0.75).*(x>=0.25);
+
+x_ref = linspace(0,1,1001);
+
+dx = 'cd'; dt = 'be';
+
+for nic = 1 : length(f_ic)
+  figure('units','normalized','outerposition',[0 0 1 1]);
+  hold on;
+  u_ref = f_ic{nic}(x_ref);
+
+  for m = 1 : length(cfl{nic})
+    [x_sol, u_sol] = adv_cauchy(N,cfl{nic}(m),c,dt,dx,f_ic{nic});
+    plot(x_sol, u_sol, cfl_spec{m}, 'LineWidth', 4, 'MarkerSize', 12);
+    legs{m} = sprintf('cfl=%f',cfl{nic}(m));
+  end
+  plot(x_ref, u_ref, 'k-', 'LineWidth', 4);
+  legs{length(cfl{nic})+1} = 'Exact';
+  xlim([0 1]);
+  ylim([-.2 1.2]);
+  xlabel('x','FontSize',36);
+  ylabel('u','FontSize',36);
+  leg = legend(legs, 'Location', 'Best');
+  set( leg, 'FontSize', 24 );
+  set( gca, 'FontSize', 24 );
+end

MATH_5343_Numerical_Solutions_to_PDE_Zeng/adv/plot_example_cd_cn.m

@@ -0,0 +1,32 @@
+c = 1.0;
+N = 80;
+cfl = {[0.5 0.99 5], [0.5 0.99 5]};
+cfl_spec = {'kx--','k*-.','ks:'};
+
+f_ic{1} = @(x)16*(x.*(1-x)).^2;
+f_ic{2} = @(x)(x<=0.75).*(x>=0.25);
+
+x_ref = linspace(0,1,1001);
+
+dx = 'cd'; dt = 'cn';
+
+for nic = 1 : length(f_ic)
+  figure('units','normalized','outerposition',[0 0 1 1]);
+  hold on;
+  u_ref = f_ic{nic}(x_ref);
+
+  for m = 1 : length(cfl{nic})
+    [x_sol, u_sol] = adv_cauchy(N,cfl{nic}(m),c,dt,dx,f_ic{nic});
+    plot(x_sol, u_sol, cfl_spec{m}, 'LineWidth', 4, 'MarkerSize', 12);
+    legs{m} = sprintf('cfl=%f',cfl{nic}(m));
+  end
+  plot(x_ref, u_ref, 'k-', 'LineWidth', 4);
+  legs{length(cfl{nic})+1} = 'Exact';
+  xlim([0 1]);
+  ylim([-.2 1.2]);
+  xlabel('x','FontSize',36);
+  ylabel('u','FontSize',36);
+  leg = legend(legs, 'Location', 'Best');
+  set( leg, 'FontSize', 24 );
+  set( gca, 'FontSize', 24 );
+end

MATH_5343_Numerical_Solutions_to_PDE_Zeng/adv/plot_example_lf.m

@@ -0,0 +1,30 @@
+c = 1.0;
+N = 80;
+cfl = {[0.99 1.08 1.15], [0.99 1.02 1.05]};
+cfl_spec = {'kx--','k*-.','ks:'};
+
+f_ic{1} = @(x)16*(x.*(1-x)).^2;
+f_ic{2} = @(x)(x<=0.75).*(x>=0.25);
+
+x_ref = linspace(0,1,1001);
+
+for nic = 1 : length(f_ic)
+  figure('units','normalized','outerposition',[0 0 1 1]);
+  hold on;
+  u_ref = f_ic{nic}(x_ref);
+
+  for m = 1 : length(cfl{nic})
+    [x_sol, u_sol] = adv_cauchy_lf(N,cfl{nic}(m),c,f_ic{nic});
+    plot(x_sol, u_sol, cfl_spec{m}, 'LineWidth', 4, 'MarkerSize', 12);
+    legs{m} = sprintf('cfl=%f',cfl{nic}(m));
+  end
+  plot(x_ref, u_ref, 'k-', 'LineWidth', 4);
+  legs{length(cfl{nic})+1} = 'Exact';
+  xlim([0 1]);
+  ylim([-.2 1.2]);
+  xlabel('x','FontSize',36);
+  ylabel('u','FontSize',36);
+  leg = legend(legs, 'Location', 'Best');
+  set( leg, 'FontSize', 24 );
+  set( gca, 'FontSize', 24 );
+end

MATH_5343_Numerical_Solutions_to_PDE_Zeng/adv/plot_example_upw_be.m

@@ -0,0 +1,32 @@
+c = 1.0;
+N = 80;
+cfl = {[0.99 2 5], [0.99 2 5]};
+cfl_spec = {'kx--','k*-.','ks:'};
+
+f_ic{1} = @(x)16*(x.*(1-x)).^2;
+f_ic{2} = @(x)(x<=0.75).*(x>=0.25);
+
+x_ref = linspace(0,1,1001);
+
+dx = 'upw'; dt = 'be';
+
+for nic = 1 : length(f_ic)
+  figure('units','normalized','outerposition',[0 0 1 1]);
+  hold on;
+  u_ref = f_ic{nic}(x_ref);
+
+  for m = 1 : length(cfl{nic})
+    [x_sol, u_sol] = adv_cauchy(N,cfl{nic}(m),c,dt,dx,f_ic{nic});
+    plot(x_sol, u_sol, cfl_spec{m}, 'LineWidth', 4, 'MarkerSize', 12);
+    legs{m} = sprintf('cfl=%f',cfl{nic}(m));
+  end
+  plot(x_ref, u_ref, 'k-', 'LineWidth', 4);
+  legs{length(cfl{nic})+1} = 'Exact';
+  xlim([0 1]);
+  ylim([-.2 1.2]);
+  xlabel('x','FontSize',36);
+  ylabel('u','FontSize',36);
+  leg = legend(legs, 'Location', 'Best');
+  set( leg, 'FontSize', 24 );
+  set( gca, 'FontSize', 24 );
+end

MATH_5343_Numerical_Solutions_to_PDE_Zeng/adv/plot_example_upw_fe.m

@@ -0,0 +1,32 @@
+c = 1.0;
+N = 80;
+cfl = {[0.99 1.08 1.1], [0.99 1.005 1.01]};
+cfl_spec = {'kx--','k*-.','ks:'};
+
+f_ic{1} = @(x)16*(x.*(1-x)).^2;
+f_ic{2} = @(x)(x<=0.75).*(x>=0.25);
+
+x_ref = linspace(0,1,1001);
+
+dx = 'upw'; dt = 'fe';
+
+for nic = 1 : length(f_ic)
+  figure('units','normalized','outerposition',[0 0 1 1]);
+  hold on;
+  u_ref = f_ic{nic}(x_ref);
+
+  for m = 1 : length(cfl{nic})
+    [x_sol, u_sol] = adv_cauchy(N,cfl{nic}(m),c,dt,dx,f_ic{nic});
+    plot(x_sol, u_sol, cfl_spec{m}, 'LineWidth', 4, 'MarkerSize', 12);
+    legs{m} = sprintf('cfl=%f',cfl{nic}(m));
+  end
+  plot(x_ref, u_ref, 'k-', 'LineWidth', 4);
+  legs{length(cfl{nic})+1} = 'Exact';
+  xlim([0 1]);
+  ylim([-.2 1.2]);
+  xlabel('x','FontSize',36);
+  ylabel('u','FontSize',36);
+  leg = legend(legs, 'Location', 'Best');
+  set( leg, 'FontSize', 24 );
+  set( gca, 'FontSize', 24 );
+end

MATH_5343_Numerical_Solutions_to_PDE_Zeng/adv/startup.m

@@ -0,0 +1 @@
+beep off;

MATH_5343_Numerical_Solutions_to_PDE_Zeng/adv2d/adv2d_cauchy.m

@@ -0,0 +1,106 @@
+function [x,y,U] = adv2d_cauchy(N1,N2,cfl,c1,c2,dt,dx)
+  % [x, y, U] = adv2d_cauchy(N1,N2,cfl,c1,c2,dt,dx)
+  % Compute the solution to the Cauchy problem:
+  %    u_t+ c1 u_x + c2 u_y=0         0 <= x, y, t <= 1
+  %    u(0,y,t) = u(1,y,t)            0 <= y, t <= 1
+  %    u(x,0,t) = u(x,1,t)            0 <= x, t <= 1
+  %    u(x,y,0) = [max(0,1-(4x-2)^2-(4y-2)^2)]^2 
+  %                                   0 <= x, y <= 1
+  % using a uniform grid with N sub-intervals.
+  % Other arguments:
+  %  dt:  'be' | 'fe'
+  %  dxx: 'upw'
+
+  x = linspace(0,1,N1+1);
+  y = linspace(0,1,N2+1);
+  h1 = 1./N1;
+  h2 = 1./N2;
+  T = 1;
+  % Determine the time step size by the CFL condition 
+  k = cfl*1/(abs(c1)/h1+abs(c2)/h2);
+
+  % Index transformation from 2D to 1D: i,j -> J
+  idx = @(i,j)((i-1)*(N2+1)+j);
+
+  % Create a one-dimensional array for the 2D initial data
+  u = zeros((N1+1)*(N2+1),1);
+  for i = 1:N1+1
+    for j = 1:N2+1
+      u(idx(i,j)) = (max(0,1-(4*x(i)-2)^2-(4*y(j)-2)^2))^2;
+    end
+  end
+
+  % Spatial discretization in x
+  A = zeros((N1+1)*(N2+1));
+  I = eye((N1+1)*(N2+1));
+  % Setup index for periodic boundary condition
+  rght = @(i)(i<N1+1)*(i+1)+(i==N1+1)*2;
+  left = @(i)(i>1)*(i-1)+(i==1)*N1;
+  abve = @(j)(j<N2+1)*(j+1)+(j==N2+1)*2;
+  belw = @(j)(j>1)*(j-1)+(j==1)*N2;
+  % - discretization of diffusion part -
+  if strcmp( dx, 'upw' ) == 1
+    if c1 > 0.0
+      for i = 1 : N1+1
+        for j = 1 : N2+1
+          J = idx(i,j); J_left = idx(left(i),j);
+          A(J,J_left) = A(J,J_left) - c1/h1;
+          A(J,J)      = A(J,J) + c1/h1;
+        end
+      end
+    else
+      for i = 1 : N1+1
+        for j = 1 : N2+1
+          J = idx(i,j); J_rght = idx(rght(i),j);
+          A(J,J)      = A(J,J) - c1/h1;
+          A(J,J_rght) = A(J,J_rght) + c1/h1;
+        end
+      end
+    end
+    if c2 > 0.0
+      for i = 1 : N1+1
+        for j = 1 : N2+1
+          J = idx(i,j); J_belw = idx(i,belw(j));
+          A(J,J_belw) = A(J,J_belw) - c2/h2;
+          A(J,J)      = A(J,J) + c2/h2;
+        end
+      end
+    else
+      for i = 1 : N1+1
+        for j = 1 : N2+1
+          J = idx(i,j); J_abve = idx(i,abve(j));
+          A(J,J)      = A(J,J) - c2/h2;
+          A(J,J_abve) = A(J,J_abve) + c2/h2;
+        end
+      end
+    end
+  else
+    disp('The diffusion part is not supported');
+  end
+  % - function handle to march in time -
+  if strcmp( dt, 'fe' ) == 1
+    mt = @(kval,uval)(I-kval*A)*uval;
+  elseif strcmp( dt, 'be' ) == 1
+    mt = @(kval,uval)(I+kval*A)\uval;
+  else
+    disp('The time-marching method is not supported');
+  end
+  % - solving the advection equation -
+  finish = 0;
+  while finish==0
+    if k>T % Check if the remaining time is smaller than k
+      k = T;
+      finish = 1;
+    end
+    u = mt(k,u); % Marching in time
+    T = T - k; % Remove time step size from remaining time.
+  end
+
+  % - rearrange the data into a 2D-array
+  U = zeros(N1+1,N2+1);
+  for i = 1:N1+1
+    for j = 1:N2+1
+      U(i,j) = u(idx(i,j));
+    end
+  end
+end