Variable step size integration

24 次查看（过去 30 天）
James Wright
James Wright 2012-5-26
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));
  2 个评论
Oleg Komarov
Oleg Komarov 2012-5-26
Supply an example of call to
ConformalSymplecticCubic2while(eps,y,tspan,x0,h)
so that we can run it and profile.
James Wright
James Wright 2012-5-26
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Oleg Komarov
Oleg Komarov 2012-5-26
You know t, htemp, and tspan(2), then you should preallocate. It will speed up things.
  2 个评论
James Wright
James Wright 2012-5-26
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
Oleg Komarov 2012-5-26
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.

请先登录，再进行评论。

John D'Errico
John D'Errico 2012-5-26
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.

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by