Richardson extrapolation for time scheme
10 views (last 30 days)
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 :
(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
Torsten
on 30 May 2022
Edited: Torsten
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)
See Also
Categories
Find more on Boundary Conditions 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!