I'm using ode45()
in matlab for some work in dynamics, using this procedure to calculate the largest Lyaponov exponent of the Lorenz system.
The process involves solving a system of differential equations starting from x0
, and comparing this with a trajectory starting very close to x0
.
At each time step, the second trajectory needs to be re-adjusted before advancing the time step, so I would like to be able to call ode45()
just once - is this possible?
A starting attempt for an alternative solution is presented here, but it doesn't work; as far as I can see, the matrices r
and R
obtained below should be similar:
% Time
ti = 0; tf = 1; res = 10;
T = linspace(ti, tf, res);
% Solve the system first time
[~,R] = ode45('tbmLorenz',T,x0);
% Alternate trajectory
r = zeros(size(R));
% Temporary cordinates
temp = zeros(3,3);
% Solve for second trajectory
for i = 2:(res-1)
% Time step
ts = T((i-1):(i+1));
% Solve three steps
[~,temp] = ode45('tbmLorenz',ts,r(i-1,:));
r_i = temp(2,:)
% Code to alter r_i goes here
% Save this
r(i,:) = r_i;
end
... but they aren't:
r =
1.0000 3.0000 4.0000
9.7011 20.6113 7.4741
29.9265 16.4290 79.0449
-5.7096 -15.2075 49.2946
-12.4917 -13.6448 44.7003
-13.6131 -13.8826 45.0346
-13.5061 -13.1897 45.4782
-13.0538 -13.0119 44.5473
-13.4463 -13.8155 44.4783
0 0 0
>> R
R =
1.0000 3.0000 4.0000
9.7011 20.6139 7.4701
29.9663 16.5049 79.1628
-5.7596 -15.2745 49.3982
-12.4738 -13.5598 44.7800
-13.5440 -13.8432 44.9084
-13.5564 -13.3049 45.4568
-13.1016 -12.9980 44.6882
-13.3746 -13.7095 44.4364
-13.7486 -13.6991 45.4092
That the last line of r
is zero is not a problem.
Any ideas? Cheers! \T