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
Ode45 is an adaptive algorithm with varying step size by definition. You can fix the step size if you want, but then you are not using ode45, you're using something else.
If you want to solve this type of problem, I suggest you use a fixed time step. This link has code for what you need.
If you insist on using ode45, then you must adjust the simulation file itself with time-interpolated inputs.