A Help understanding this extrapolation technique

If you Taylor expand out the function in time, you get for time step h:
##r(t+h) = r(t) + h r'(t) + \frac{h^2}{2} r''(t) + \frac{h^3}{3!} r'''(t)...\\
r(t - h) = r(t) - h r'(t) + \frac{h^2}{2} r''(t) -\frac{h^3}{3!} r'''(t)... \\
r(t-2h) = r(t) - 2h r'(t) + \frac{4h^2}{2} r''(t) -\frac{8h^3}{3!} r'''(t)... ##
##3r(t) - 3r(t-h) + r(t-2h) = r(t) +h r'(t) + \frac{h^2}{2} r''(t) -\frac{5h^3}{6} r'''(t)... ##
So you can see this scheme is accurate to ##\mathcal{O}(h^3). ##

The benefit of this method is that it looks like an semi-implicit method, which should become stable regardless the sample rate. The accuracy will still be dependent on sample rate, but the stability is much improved compared to forward difference methods.

 

In the attachment, the equation marked (6.46) says
##\frac{dr_1}{dt} = f(\hat r _1, t_1).##
Do you know that function? If not, you can approximate the derivatives simply, but you lose some accuracy.


Here is a rough code for the approximation:
-------------------------------------------------------------

function R = bubble(r0,r1,r2, N, dt)
% the output of will be a Nx2 matrix with the first column as the radius
% and the second column the first derivative in time. Inputs are the radius
% at three known times (r1,r2,r3), the number of timesteps (N),
% the step size (dt), and the function f.

%initialize with initial conditions
rnow = r2;
rnowless1 = r1;
rnowless2 = r0;

R = zeros(N+3,1);
R(1) = r0;
R(2) = r1;
R(3) = r2;
%timestep routine
for step = 1:N
% guess at rnext.
rhat = 3*rnow -3*rnowless1 + rnowless2; %eq 6.46 from text
rbardot = (rhat-rnowless1)/(2*dt); %this is the same as [(rhat-rnow)/dt + (rnow-rnowless1)/dt]/2
rnext = rnow + rbardot*dt;
R(step+3) = rnext; %save radius data.
% reset for new step
rnowless2 = rnowless1;
rnowless1 = rnow;
rnow = rnext;
end

plot(0:dt:(N+2)*dt,R);
end

 

That seems to work, remember that you are using the approximation ##\hat r_1 ## to approximate the derivative, so plug that estimate (based on prior r data) into the formula for dr/dt.
Once you have that estimate for dr/dt, you apply the averaging and approximation in 6.48 and 6.49.

Assuming you know the previous time data ##r_{i-1},r_{i-2},r_{i-3}##:
1. Set dt timestep.
2. Set pressure rate.
3. Calculate approximation ##\hat r_i## based on previous r information.
4. Use formula with pressure and ##\hat r_i## or ## r_{i-1}## to get ##\frac{dr_i}{dt}## and ##\frac{dr_{i-1}}{dt}##
5. Use equation 6.48 to get average rate of change for the timestep.
6. Use equation 6.49 to approximate ##r_i##.

 

Based on the example, it looks like you would want to have the radius at the first 3 times. You can surely assume that r = 0 or r = r0 for times prior to start.
Yes, I can look at a script for you.

 

If you know ##r_{i-1}##, you just plug it in. There is nothing to solve.
For the current (i) timestep, you don't know ##r_i##, so you use the ##\hat r_i## that you already computed.