Richardson extrapolation for time scheme
Show older comments
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 :
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.
14 Comments
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
on 27 May 2022
Here is a survey article on 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
on 28 May 2022
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
on 28 May 2022
Karl Zero
on 29 May 2022
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.
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
on 29 May 2022
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
on 30 May 2022
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
Answers (0)
Categories
Find more on Linear Algebra in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
