Variable step size integration

Hi, I've written a multi Symplectic integrator and want to make the program vary it's step size depending on the dynamics. I have written a program which I believe does this, but it takes much longer than the original, and the longer the time I want to integrate over, the worse it gets. i.e. time 0:50 it takes 33 times as long, 0:500 it took 147 times as long, 0:1000 I gave up waiting after 10 minutes. This is a big problem as I wish to be able to integrate from 0:500000 and longer. If anyone can see how I can either speed my code up, or suggest a different, faster method, it'd be a great help. See below for both codes. Thanks,

First here is my initial integrator

function [U] = ConformalSymplecticCubic2(eps,y,tspan,x0,h)
% eps is a parameter (0.1 in my case)
% y is another parameter which represents dissipation
%   this should be taken to be small, i.e 0.0005 or less
% tspan is a vector [tstart,tfinal]
% x0 is a vector for the initial conditions [x,x']
% h is the step size to be used, I use 0.01
tic
N = (tspan(2)-tspan(1))/h;
%setting vector sizes
U = zeros(N,2);
t = zeros(N,1);
t(1) = tspan(1);
%assigning inital values
U(1,1) = x0(1);
U(1,2) = x0(2);
%Integration
for i=2:N
    U(i,2) = exp(-h*y)*U(i-1,2) - h*(1+eps*cos(t(i-1)))*(U(i-1,1))^3;
    U(i,1) = U(i-1,1) + (h)*U(i,2);
    t(i) = t(i-1) + h;
end
toc
plot(U(:,1),U(:,2));

And here is the one with the variable step size

function [U] = ConformalSymplecticCubic2while(eps,y,tspan,x0,h)
tic
Tol = 10^(-5);
cmax = 20;
Utemp = zeros(2,2);
Utemp2 = zeros(3,2);
ttemp = zeros(3,1);
t = tspan(1);
U(1,1) = x0(1);
U(1,2) = x0(2);
D = exp(-h*y);
i = 2;
while t < tspan(2)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Integrate for original step size %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    U(i,2) = D*U(i-1,2) - h*(1+eps*cos(t))*(U(i-1,1))^3;
    U(i,1) = U(i-1,1) + (h)*U(i,2);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Finished initial integration now %
    %  checking for smaller step size  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    c = 0;
    dist = 1;
    htemp2 = h;
    Utemp(1,1) = U(i-1,1);
    Utemp(1,2) = U(i-1,2);
    Utemp2(1,1) = Utemp(1,1);
    Utemp2(1,2) = Utemp(1,2);
    Utemp2(2,1) = U(i,1);
    Utemp2(2,2) = U(i,2);
    ttemp(1) = t;
    while dist > Tol && c <= cmax
        %Makes distance of next integration
        %half of this time round (should it
        %be needed)
        Utemp(2,1) = Utemp2(2,1);
        Utemp(2,2) = Utemp2(2,2);
        c = c+1;
        htemp = htemp2;
        %setting up initial temporary constants
        norm1 = sqrt(Utemp(2,1)^2 + Utemp(2,2)^2);
        htemp2 = htemp/2;
        Dtemp = exp(-htemp2*y);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Integrating for half step size of %
        %        previous integration       %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        for k = 2:3
            Utemp2(k,2) = Dtemp*Utemp2(k-1,2) - htemp2*(1+eps*cos(ttemp(k-1)))*(Utemp2(k-1,1)^3);
            Utemp2(k,1) = Utemp2(k-1,1) + htemp2*Utemp2(k,2);
            ttemp(k) = ttemp(k-1) + htemp2;
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Checking distance between last 2 %
        %           Integrations           %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        norm2 = sqrt(Utemp2(3,1)^2 + Utemp2(3,2)^2);
        dist = abs(norm1 - norm2);
        if c == cmax
            S = 1;
        end
    end
    U(i,1) = Utemp(2,1);
    U(i,2) = Utemp(2,2);
    t = t+htemp;
    i = i +1;
end
toc
plot(U(:,1),U(:,2));

[Oleg Komarov]

Supply an example of call to

ConformalSymplecticCubic2while(eps,y,tspan,x0,h)

so that we can run it and profile.

[James Wright]

ConformalSymplecticCubic2while(0.1,0.0005,[0,500],[0.535,0.791],0.01);

Similarly

ConformalSymplecticCubic2(0.1,0.0005,[0,500],[0.535,0.791],0.01);

will run my previous program, for comparison of speed.

Answer 1

[Oleg Komarov]

You know t, htemp, and tspan(2), then you should preallocate. It will speed up things.

[James Wright]

but htemp changes, so I don't know how many time steps will be taken, so I can't preallocate. I have preallocated where I know sizes

[Oleg Komarov]

Not trivial, but you could preallocate for known h, then whenever htemp2 halves add additional space (that can be calculated since htemp2 is a known function of h). It will be less expensive to dynamically preallocate than to grow U.

Now, since I stated "not trivial" and not being an expert in numerical integration I would recommend looking around for known algorithms on how to approach this problem before stepping into dynamic preallocation.

Answer 2

[John D'Errico]

It looks like you have not learned why it is good to preallocate some arrays. Instead, they keep growing. Perhaps this is time to learn that lesson.

[James Wright]

please see above