Richardson extrapolation for time scheme

Question

Karl Zero 2022-5-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1728420-richardson-extrapolation-for-time-scheme

编辑： Torsten 2022-5-30

在 MATLAB Online 中打开

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.

14 个评论
显示 12更早的评论隐藏 12更早的评论

Karl Zero 2022-5-28

在 MATLAB Online 中打开

@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

Karl Zero 2022-5-29

编辑：Karl Zero 2022-5-29

在 MATLAB Online 中打开

@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 2022-5-29

编辑：Torsten 2022-5-29

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.

Torsten 2022-5-29

编辑：Torsten 2022-5-30

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 2022-5-30

在 MATLAB Online 中打开

@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 2022-5-30

编辑：Torsten 2022-5-30

在 MATLAB Online 中打开

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