Finite Difference Centre Time Centre Space

Hi, I'm trying to integrate the 1D linear advection equation using a Centred-in-time, Centred-in-space method. I was given some code by my professor for a FTCS scheme (which worked) and I've tried to adapt it using the CTCS method but the problem is that phim(j) which is used to calculate phip(j) is not defined anywhere in the code and I don't know how/ where to define it.

I will attach the code I have so far:

function val=icA(x)
%Initial Condition for Advection Problem
if (0.25 <= x) & (x <= 0.75)
    val = 0.5*(1+cos(4*pi*(x-0.5)));
else
    val = 0;
end;

This describes the initial conditions for the function

%Set up plots to solve linear advection equation using a CTCS scheme
u = 2.0; %Set constant velocity value
nx=51; %Number of Grid Points
dx = 1.0/(nx-1); %Creates grid length of 0.02 (1/50)
x = (0:1:nx-1).*dx; %Discretise in space for 50 equally-spaced grid intervals
phi=zeros(nx,1); %Create vector that stores values of phi for each grid point with nx rows, 1 column
phip=zeros(nx,1); %Create vector that stores integrated values of phi with nx rows, 1 column
dt = 0.005; %Create timestep of 0.005
nu = (u*dt)/dx;
disp(nu);
%Set Initial Condition
for j = 1:nx
    phi(j)=icA(x(j));
end;
plot(x,phi); %Plot x against phi
phi_i_c=phi;
nstep=101;
%Loop over steps
for istep=1:nstep
    %Loop over grid points
    for j=1:nx
        jm = j-1; %Initial evaluation of jm = j-1
        if jm < 1
        jm = jm + nx - 1; %Modify jm to account for periodic BCs.
    end;
    jp = j+1; %Initial evaluation of jp = j+1
    if jp > nx
        jp = jp - nx + 1; %Modify jp to account for periodic BCs
    end; 
    %phim(j) = ????
    phip(j)=phim(j)-nu*(phi(jp)-phi(jm))

end; %j
phi=phip;
end; %istep
%
hold on
plot(x,phi,'r')

This outlines the CTCS scheme I have curated thus far.

Essentially I am just wondering how to embed phim(j) into the code. Any help in the right direction would be appreciated.

Answer 1

[VBBV]

%Set up plots to solve linear advection equation using a CTCS scheme
u = 2.0; %Set constant velocity value
nx=51; %Number of Grid Points
dx = 1.0/(nx-1); %Creates grid length of 0.02 (1/50)
x = (0:1:nx-1).*dx; %Discretise in space for 50 equally-spaced grid intervals
phi=zeros(nx,1); %Create vector that stores values of phi for each grid point with nx rows, 1 column
phip=zeros(nx,1); %Create vector that stores integrated values of phi with nx rows, 1 column
phim=zeros(nx,1);
dt = 0.005; %Create timestep of 0.005
nu = (u*dt)/dx;
disp(nu);

    0.5000

%Set Initial Condition
for j = 1:nx
    phi(j)=icA(x(j));
end;
subplot(211)
plot(x,phi); %Plot x against phi
phi_i_c=phi;
nstep=101;
%Loop over steps
for istep=1:nstep
    %Loop over grid points
    for j=1:nx
        jm = j-1; %Initial evaluation of jm = j-1
        if jm < 1
        jm = jm + nx - 1; %Modify jm to account for periodic BCs.
    end;
    jp = j+1; %Initial evaluation of jp = j+1
    if jp > nx
        jp = jp - nx + 1; %Modify jp to account for periodic BCs
    end; 
    %phim(j) = ????
    phip(j)=phim(j)-nu*(phi(jp)-phi(jm));

end; %j
phi(:,istep)=phip;  % use phip here 
end; %istep
%
subplot(212)
plot(x,phi,'r')

function val=icA(x)
%Initial Condition for Advection Problem
if (0.25 <= x) & (x <= 0.75)
    val = 0.5*(1+cos(4*pi*(x-0.5)));
else
    val = 0;
end;
end

Define the phim array at beginnig like other arrays and compute for every step

[Angus K]

Thanks you're correct that I missed defining the phim array at the start but how can phip(j) be evaluated when phim(j) is not defined anywhere in the script?

[VBBV]

phim requires preallocation of array at beginning

phim=zeros(nx,1);   % allocate zeros array at beginning 

Then compute phip as below

 phip(j)=phim(j)-nu*(phi(jp)-phi(jm));  % this line 
 phim(j) = %use the equation 

is phim changing for every step ? is there any equation for phim ?

[Angus K]

Well, phim as I intend to define it is just the value of phi two time iterates before phip if that helps. It just follows from the centre time centre space (CTCS) integration scheme if you’re familiar. I’ve attached an image if not. You essentially just rearrange for Phi^(n+1) and that’s your scheme for integration

Answer 2

[Torsten]

Never ever use centered-in-space for the advection equation !

I commented out the settings for the centered-in-space scheme (nu and phip) and included the stable upwind scheme.

Greater modifications will be necessary if you want to include "centered-in-time". In the code below, "backward in time" is used.

%This describes the initial conditions for the function
%Set up plots to solve linear advection equation using a CTCS scheme
u = 2.0; %Set constant velocity value
nx=51; %Number of Grid Points
dx = 1.0/(nx-1); %Creates grid length of 0.02 (1/50)
x = (0:1:nx-1).*dx; %Discretise in space for 50 equally-spaced grid intervals
phi=zeros(nx,1); %Create vector that stores values of phi for each grid point with nx rows, 1 column
phip=zeros(nx,1); %Create vector that stores integrated values of phi with nx rows, 1 column
dt = 0.005; %Create timestep of 0.005
%nu = (u*dt)/(2*dx);   % nu for centered-in-space
nu = (u*dt)/dx;        % nu for upwind
disp(nu);
%Set Initial Condition
for j = 1:nx
    phi(j)=icA(x(j));
end;
plot(x,phi); %Plot x against phi
phi_i_c=phi;
nstep=51;
%Loop over steps
for istep=1:nstep
    %Loop over grid points
    for j=1:nx
        jm = j-1; %Initial evaluation of jm = j-1
        if jm < 1
            jm = jm + nx - 1; %Modify jm to account for periodic BCs.
        end;
        jp = j+1; %Initial evaluation of jp = j+1
        if jp > nx
            jp = jp - nx + 1; %Modify jp to account for periodic BCs
        end; 
        %phip(j)=phi(j)-nu*(phi(jp)-phi(jm))   %phip for centered-in-space
        phip(j) = phi(j)-nu*(phi(j)-phi(jm));  %phip for upwind
    end; %j
    phi=phip;
    hold on
    plot(x,phi,'r')
end; %istep
%
%hold on
%plot(x,phi,'r')

function val=icA(x)
%Initial Condition for Advection Problem
if (0.25 <= x) & (x <= 0.75)
    val = 0.5*(1+cos(4*pi*(x-0.5)));
else
    val = 0;
end;
end

[Angus K]

Hey Torsten.

Thanks for this. Problem is I’ve also never used CTCS before but my professor specifically asks for the question to be solved like this..

[Torsten]

You will have to ask your professor with which scheme you should do the first integration step (from time t=0 to time t = dt) since the suggested scheme cannot be applied for this starting step.

It's the same as with a recursion

x_(n+2) = x_n + x_(n+1).

You must know x_0 and x_1 before you can start.

[Angus K]

Hi Torsten,

Just had a read through the lecture notes. I did think this too and it suggests using FTCS to get started.