Richardson extrapolation for time scheme

I am solving an equation by using finite volume method :

I am writing a linear system in order to solve it by using a forward euler scheme in time. The scheme is written as :

It is a semi-implicit scheme as the derivative of the Q term is writing in the time (n+1). We end up with a linear system to solve as :

My question is now : how to write the code that is related to this scheme in order to use a Richardson extrapolation ?

I have read some papers in order to understand the method but I have some misunderstandings... For instance, if we evaluate the value of the function u at time step h and the value of this same function in the time step 2h we can write the exact value of that one is trying to compute as :

(k is the order of approximation of the method)

I don't know how I can implement it in my code since I am basically computing a time loop solving the linear system at each time step. Moreover : the exact solution of my problem is unknown...

Nx = 100;                    % Number of grid points
L = 10;                      % Length of interval
x = linspace(0, L, Nx+1);    % Grid points
dx = x(2) - x(1);            % Spatial step
t = 0;                       % Initial time
dt = 0.1;                    % Time step
t_final = 30;                % Final time
Nt = t_final/dt;             
h_new   = zeros(1,Nx+1);     % h at new time
h_old = zeros(1,Nx+1);       % h at old time
c =  dt/(dx^2);              % Constant
epsilon = 1e-7;              % condition

% Initialization of the linear system
A = zeros(Nx+1, Nx+1);
b = zeros(Nx+1,1);
solution = zeros(Nt,Nx+1);

% Stationnary condition
for i=1:Nx+1
    h_old(i)= (1-((x(i)/L)).^2).^2;  
end

%% Solving AU=b 

% Time loop
while t<=t_final
    % Matrix elements
    for i=2:Nx

        A(i,i-1) = -c*((h_old(i)+h_old(i-1))/2)^3;
        A(i,i+1) = -c*((h_old(i)+h_old(i+1))/2)^3;
        A(i,i)   = 1 +c*((h_old(i)+h_old(i+1))/2)^3+c*((h_old(i)+h_old(i-1))/2)^3;

    end

    A(1,1) = (1+2*c*((h_old(1)+h_old(2))/2)^3);
    A(1,2) = -2*c*((h_old(1)+h_old(2))/2)^3;
    A(Nx+1,Nx+1) = 1;


    % b vector elements
    for i=1:Nx
        b(i) = h_old(i);
    end
    b(Nx+1) = epsilon;
    h_new(:) = linsolve(A, b);

    % Update of h_old.
    h_old(:) = h_new;

    % Time iteration.
    t = t+dt;
end

Could someone help me to do so please ? And understand the implementation of the method.

[Torsten]

If you know the order of your method above, you'll have to solve the problem several times with Nx = 100, Nx = 200, Nx = 400, e.g., and extrapolate the results to dx=0 with the formula you've written.

Are you sure about the boundary condition at x=L ? The scheme produces strange curves for greater values of tfinal.

[Karl Zero]

@Torsten I don't understand your point : you mean I have to compute the problem for multiple number of points ? But how do i make it through a loop ? For instance calculating for dt than for 2*dt. The boundary condition at x=L is defined on the function itself so I think that is the correct way to do it. But my main concern here is to learn the Richardson extrapolation to be honest.

[Torsten]

Here is a survey article on Richardson extrapolation:

https://en.wikipedia.org/wiki/Richardson_extrapolation

If you have a specific question, you may come back here.

The loop in which you can calculate a solution of your problem for different numbers of grid points (and thus dx values) can start just at the beginning of your code and end at the end of your code.

In the present form, "solution" is a matrix of size (Nt,Nx+1).

In the case you work with different grids, you should introduce a cell array SOL in which you save "solution" for different values for Nx:

Nx = 100;
for igrid = 1:ngrid
  % your code
  Nx = 2*Nx;
  SOL{igrid} = solution;
end

[Karl Zero]

@Torsten I don't understand why you advise to compute for multiple values of Nx since I want the Richardson method only for the time integration ? For example : the formula written in my post is for 2 time steps dt and 2*dt.

Thank you for the article provided : but I ended up with the same misunderstanding... For example : they provide a code to solve an ODE and they choose to calculate on 2 time steps : dt and dt/2 and they construct a loop to compute the matrix elements :

    for j = 1 : i     
        
        % Go across this current (i+1)-th row until the diagonal is reached
        % To compute A(i + 1, j + 1), which is the next Richardson extrapolate, 
        % use the most recently computed value (i.e. A(i + 1, j)) and the value from the
        % row above it (i.e. A(i, j)).
     
        A(i + 1, j + 1) = ((4^j).*A(i + 1, j) - A(i, j))/(4^j - 1);
    end

For what I have understood : they compute the first element of the matrix using a trapezoidal method :

A(1, 1) = Trapezoidal(f, tStart, tEnd, h, y0)

Then they compute the second row by using the second time step (dt/2) :

for i = 1 : maxRows - 1 % Starting at i = 1, iterate at most maxRows - 1 times
    h = h/2             % Halve the previous value of h since this is the start of a new row.
    
    % Starting filling row i+1 from the left by calling the Trapezoidal function with this new smaller step size
    A(i + 1, 1) = Trapezoidal(f, tStart, tEnd, h, y0)

I am not sure about how this can be applied on my own problem since my first row elements depend on time step by the parameter

c = dt/(dx^3);

% First row elements

    A(1,1) = (1+2*c*((h_old(1)+h_old(2))/2)^3);
    A(1,2) = -2*c*((h_old(1)+h_old(2))/2)^3;

Then the computing of the matrix elements would also depend on "c" :

for i=2:Nx
    
    A(i,i-1) = -c*((h_old(i)+h_old(i-1))/2)^3;
    A(i,i+1) = -c*((h_old(i)+h_old(i+1))/2)^3;
    A(i,i)   = 1 +c*((h_old(i)+h_old(i+1))/2)^3+c*((h_old(i)+h_old(i-1))/2)^3;

end

Does that mean that I am suppose to write something like :

Nx = 100; 
dt = 0.1;
tolerance = 10^-11

A = zeros(Nx+1, Nx+1);

step = @(k) k/(dx^2);
c = step(dt);

A(1,1) = (1+2*c*((h_old(1)+h_old(2))/2)^3);
A(1,2) = -2*c*((h_old(1)+h_old(2))/2)^3;


for i = 2 : Nx  
    dt = dt/2
    c = step(dt);
    
    
    A(i,i-1) = -c*((h_old(i)+h_old(i-1))/2)^3;
    A(i,i+1) = -c*((h_old(i)+h_old(i+1))/2)^3;
    A(i,i)   = 1 +c*((h_old(i)+h_old(i+1))/2)^3+c*((h_old(i)+h_old(i-1))/2)^3;
    
    for j = 2 : Nx       
        A(i + 1, j + 1) = ((4^j).*A(i + 1, j) - A(i, j))/(4^j - 1);
    end
    
    if (absoluteValue(A(i + 1, i + 1) - A(i, i)) < tolerance)  

    break                                                  
    end
end

I dont know if I am clear in my formulation sorry

[Torsten]

I don't understand why you advise to compute for multiple values of Nx since I want the Richardson method only for the time integration ?

Richardson extrapolation for time integration only makes sense if it is done together with an adaptive stepsize control in time. Making a time step dt and two time steps dt/2 give you two approximations for your unknown function in t+dt. The difference of these two approximations give you information about the local error of your time integration. If the error is above a prescribed error level, you reduce the stepsize, if it is below a prescribed error level, you enlarge the stepsize. Coupling this stepsize control with Richardson extrapolation results in a time integration method that is one order better in time than the one you use.

But in the realm of partial differential equations, I only know of the application of Richardson extrapolation for the spatial discretization error, not the time integration error. And that's where calculations on different grids come into play.

[Karl Zero]

@Torsten I think that this is what they try to do in the article you provided : don't you think that my implementation is correct ? I am not sure about the "update" of the "c" parameter of my code and also the way that my loop is constructed...

[Karl Zero]

@Torsten Do you know why the Richardson update is written as :

for j = 2 : Nx       
    A(i + 1, j + 1) = ((4^j).*A(i + 1, j) - A(i, j))/(4^j - 1);
end

I don't get why is it 4^j ? As the power should be fixed by the order of magnitude of the method ?

[Torsten]

Your matrix A in which you set up the linear system

A*h = b

to be solved in each time step to get h at t=t+dt has absolutely nothing to do with the matrix A in the article.

A in the article is the Richardson extrapolation matrix.

Your code would be the function called "Trapezoidal(f, tStart, tEnd, dt, y0)" which integrates the system in question from tStart to tEnd with a stepwise dt and initial condition y0. The function that you had to call at this place could be called "Euler_implicit(f, tStart, tEnd, dt, y0)". But note that your code does not return a single scalar as the code in the Wikipedia article, but a vector of size Nx. Thus the Richardson matrix you had to use would be of size A(Nx+1,MaxRows,MaxRows).

The trapezoidal rule has a development of the local discretization error in dt^2. That's where the 4^j instead of the 2^j in the update stems from.

If you want an adaptive time integrator for your PDE that tries to integrate within a prescribed error tolerance, use the method-of-lines and MATLAB's ode15s.

[Karl Zero]

@Torsten Do you think that I am building the function in the correct way ? I am not sure about several parameters as compared to the article. The function f should be the solved function so I think that it is the h_old of my script ? And since I don't have a value for the initial position : can I choose to not specify it ? Also do I have to specify my matrix A as an output of the function ?

I am mainly trying to procede as the article : it allows me to better understand this scheme.

Thank you a lot for your help.

Nx = 100;                    % Number of grid points
L = 10;                      % Length of interval
x = linspace(0, L, Nx+1);    % Grid points
dx = x(2) - x(1);            % Spatial step
t = 0;                       % Initial time
dt = 0.1;                    % Time step
t_final = 30;                % Final time
Nt = t_final/dt;             
h_new   = zeros(1,Nx+1);     % h at new time
h_old = zeros(1,Nx+1);       % h at old time
epsilon = 1e-7;              % condition

% Stationnary condition
for i=1:Nx+1
    h_old(i)= (1-((x(i)/L)).^2).^2;  
end

%% Solving AU=b 
function h_new = Euler(h_old,t,t_final,dt)

% Initialization of the linear system
A = zeros(Nx+1, Nx+1);
b = zeros(Nx+1,1);

c =  dt/(dx^2);

% Time loop
while t<=t_final
    % Matrix elements
    for i=2:Nx
        A(i,i-1) = -c*((h_old(i)+h_old(i-1))/2)^3;
        A(i,i+1) = -c*((h_old(i)+h_old(i+1))/2)^3;
        A(i,i)   = 1 +c*((h_old(i)+h_old(i+1))/2)^3+c*((h_old(i)+h_old(i-1))/2)^3;
    end
    A(1,1) = (1+2*c*((h_old(1)+h_old(2))/2)^3);
    A(1,2) = -2*c*((h_old(1)+h_old(2))/2)^3;
    A(Nx+1,Nx+1) = 1;
    % b vector elements
    for i=1:Nx
        b(i) = h_old(i);
    end
    b(Nx+1) = epsilon;
    h_new(:) = linsolve(A, b);
    % Update of h_old.
    h_old(:) = h_new;
    % Time iteration.
    t = t+dt;
end
end

[Torsten]

You start at t = 0 and choose an initial stepsize dt.

From the function of the article, you call your Euler-function with step sizes dt, dt/2, dt/4,... until you achieve convergence, i.e. norm(A(:,i+1,i+1)-A(:,i,i)) < tolerance where ":" is the solution of your PDE in all Nx+1 grid points.

Here, usually the step size control comes into play, i.e. for the next step in time, you would change the dt depending on whether the Richardson extrapolation needed many or few steps to converge. I think this will be too difficult for you to implement.

So for the next step, you continue to use dt as step size in time, set tStart = dt and tEnd = 2*dt, call your Euler-function from the function in the article with h_old = A(:,i+1,i+1) and so on until you reach t_final = 30.

"A" as output of your funtion is not needed - you only need h_new calculated in t=tEnd to the step size 2^(-i)*dt.

What you mean by "I don't have a value for the initial position" I don't know. Your initial values are always tStart as the tEnd of the last step and h_old = A(:,i+1,i+1) as the solution for h from the last step.

[Karl Zero]

@Torsten

So for the next step, you continue to use dt as step size in time, set tStart = dt and tEnd = 2*dt, call your Euler-function from the function in the article with h_old = A(:,i+1,i+1) and so on until you reach t_final = 30.

• Since my output is a vector of size Nx : I don't understand what should be the first element of the A matrix of the article in my case :

%Compute the top left element of the matrix. The first row of this (lower triangular) matrix has now been filled.
A(1, 1) = Trapezoidal(f, tStart, tEnd, h, y0)

If I have correctly understood it should be :

%Compute the top left element of the matrix. The first row of this (lower triangular) matrix has now been filled.
A(:,1, 1) = Trapezoidal(f, tStart, tEnd, h, y0)

Correct ?

• Also : why do you suggest to set tStart = dt and tEnd=2*dt ?It is not what it is done on the article : I don't understand this part of your reasoning.

Best regards @Torsten

[Torsten]

Correct ?

Correct. A should be of size (Nx+1,maxRows,maxRows). I guess your output vector is of size (Nx+1), not Nx.

Also : why do you suggest to set tStart = dt and tEnd=2*dt ?It is not what it is done on the article : I don't understand this part of your reasoning.

As I already said: Richardson extrapolation is always used in combination with a stepsize control. So you start with a stepsize dt1. You apply Richardson extrapolation with stepsizes dt1/2, dt1/4 etc. until you achieve convergence. At this point, you've reached t = dt1. Now you compute the stepsize dt2 for the next step such that a prescribed error tolerance is met. You again apply Richardson extrapolation with stepsizes dt2/2, dt2/4,... until convergence is achieved. At this point, you've reached t=dt1+dt2. This process is continued until you've reached tFinal.

Since you will most probably not apply stepsize control, one can come with the idea to integrate as it is done in the article: just start with a stepsize dt = Tfinal = 30. But that's nonsense since the method will immediately diverge. So my advice: Choose dt not too large for that Richardson extrapolation still converges. Thus the extrapolation loop in the article code will be run with an additional outer loop in which you restart extrapolation after each successful step from (i-1)*dt to i*dt. To determine a reasonable dt, I suggest you run your original code with several time step sizes and see up to which stepsize you still get reasonable profiles for h over time. The maximum dt for which this is the case could be a good choice for the dt chosen in your extrapolation code.

By the way: You use a mixture of implicit and explicit Euler as your time integration scheme. Do you even know (for that the application of Richardson extrapolation is correct and justified) what the order of this method is and if there exists the necessary development of the local discretization error ?

[Karl Zero]

@Torsten it may seem dumb but in fact the result plotted by the main code is a matrix nammed solution which is :

solution = zeros(Nt,Nx+1 ); 

because there is a solution calcuted on Nx+1 points for each time step (in the while loop )

So for instance the line :

M(:,1, 1) = Euler_implicit(t, t_final, dt ); 

returns an error of size (which is legit ) :

Unable to perform assignment because the size of the left side is 101-by-1 and the size of the right side is 300-by-101. 

However, i don't understand how to avoid it ....

[Torsten]

You only need h for t = t_final in the Richardson function. Thus "solution" should be initialized as zeros(Nx+1,1) in Euler_implicit.

Maybe of interest for you:

format long
tStart = 0;          % Starting time
tEnd = 5;            % Ending time
f = @(t,y) -y^2;     % The derivative of y, so y' = f(t, y(t)) = -y^2
% The solution to this ODE is y = 1/(1 + t)
y0 = 1;              % The initial position (i.e. y0 = y(tStart) = y(0) = 1)
tolerance = 10^-11;  % 10 digit accuracy is desired

maxRows = 20;                % Don't allow the iteration to continue indefinitely
initialH = tEnd - tStart;    % Pick an initial step size
haveWeFoundSolution = false; % Were we able to find the solution to within the desired tolerance? not yet.

h = initialH;

% Create a 2D matrix of size maxRows by maxRows to hold the Richardson extrapolates
% Note that this will be a lower triangular matrix and that at most two rows are actually
% needed at any time in the computation.
A = zeros(maxRows, maxRows);

%Compute the top left element of the matrix. The first row of this (lower triangular) matrix has now been filled.
A(1, 1) = Trapezoidal(f, tStart, tEnd, h, y0)

% Each row of the matrix requires one call to Trapezoidal
% This loops starts by filling the second row of the matrix, since the first row was computed above
for i = 1 : maxRows - 1  % Starting at i = 1, iterate at most maxRows - 1 times
    h = h/2;             % Halve the previous value of h since this is the start of a new row.
    
    % Starting filling row i+1 from the left by calling the Trapezoidal function with this new smaller step size
    A(i + 1, 1) = Trapezoidal(f, tStart, tEnd, h, y0);
    
    for j = 1 : i     % Go across this current (i+1)-th row until the diagonal is reached
        % To compute A(i + 1, j + 1), which is the next Richardson extrapolate, 
        % use the most recently computed value (i.e. A(i + 1, j)) and the value from the
        % row above it (i.e. A(i, j)).
        
        A(i + 1, j + 1) = ((4^j).*A(i + 1, j) - A(i, j))/(4^j - 1);
    end
    
    % After leaving the above inner loop, the diagonal element of row i + 1 has been computed
    % This diagonal element is the latest Richardson extrapolate to be computed.
    % The difference between this extrapolate and the last extrapolate of row i is a good
    % indication of the error.
    if (abs(A(i + 1, i + 1) - A(i, i)) < tolerance)  % If the result is within tolerance
        %print('y = ', A(i + 1, i + 1))           % Display the result of the Richardson extrapolation
        A(i+1,i+1)
        haveWeFoundSolution = true;
        break                                                  % Done, so leave the loop
    end
end

if (haveWeFoundSolution == false)   % If we were not able to find a solution to within the desired tolerance
    print("Warning: Not able to find solution to within the desired tolerance of ", tolerance);
    print("The last computed extrapolate was ", A(maxRows, maxRows))
end


function y = Trapezoidal (f,tStart,tEnd,h,y0)
  t = tStart:h:tEnd;
  n = numel(t)
  y = y0;
  yold = y;
  for i = 2:n
    fun = @(y)y - yold - 0.5*h*(f(t(i-1),yold)+f(t(i),y));
    y = fsolve(fun,yold);
    yold = y;
  end
end