Use ODE45 solution as new input to the main code

I'm trying to calculate the tractrix curve of one object being pulled by another.

Suppose the first object moves along an arc of a circle centered at the origin and radius 5:

function [X, Xs, Y, Ys] = child(t);
%Child
X  =  5*cos(t);  Y  =  5*sin(t);
Xs = -5*sin(t);  Ys =  5*cos(t);

Where t=0:0.01:pi/2

The starting position of the second object is X = 5, Y = -3.

To calculate the tractrix curve I use this code found on the net:

% main1.m 
clear all; close all; clc;
y0 = [5 -3]';
t = 0:0.01:pi/2;
[t y] = ode45('f',t,y0);


clf; hold on;
axis([-6 6 -6 6]);
axis('square');
plot(y(:,1),y(:,2));


[X, Xs, Y, Ys] = child(t);
plot(X,Y,':')
hold off;

Where f function is define by:

function zs = f(t,z)
%
[X Xs Y Ys] = child(t);
v =[Xs; Ys];
w =[X-z(1); Y-z(2)];
w = w/norm(w);
zs = (v'*w)*w;

The result is show in figure 1:

What I want to do now is to use the solution to calculate the tractrix curve of a third object connected to the second with initial position X = 5, Y = -6. How can I do that? How can I use the ODE45 solution as new input?

Answer 1

[Cris LaPierre]

Perhaps I don't understand the problem, but it seems like all you need to solve for the position of the third object is t and its initial position.

y02 = [5; -6];
[t2 y2] = ode45('f',t,y02);

[Andrea Scodro]

Imagine you have the situation show in the picture:

The path of the 1st hook is an arc of a circle of radius 5, I want to calculate the path of the other 2 hooks. The code I posted should return the path of the 2nd hook. For the path of the 3rd hook I have to iterate the code using 2nd hook path as input.

The initial positions:

1st hook : (5,0)

2nd hook : (5,-3)

3rd hook: (5,-6)

PS: Sorry for my bad English

[Cris LaPierre]

You English is great. My point is that the inputs to ode45 are a function, which is the same for all 3 objects, a time interval or span, which is also the same for all 3 objects, and the initial postion of each object, which is unique for each object. You do not need the results from one object to solve for the position of another.

t = 0:0.01:pi/2;

% Oject 1
[X, Xs, Y, Ys] = child(t);
plot(X,Y,':','DisplayName','Obj1')

% Object 2
y01 = [5; -3];
[t y] = ode45(@f,t,y01);

% Object 3
y02 = [5; -6];
[t y2] = ode45(@f,t,y02);

hold on;
plot(y(:,1),y(:,2),'DisplayName','Obj2');
plot(y2(:,1),y2(:,2),'DisplayName','Obj3');
hold off;
axis([-6 6 -6 6]);
axis('square');
legend

function zs = f(t,z)
%
[X Xs Y Ys] = child(t);
v =[Xs; Ys];
w =[X-z(1); Y-z(2)];
w = w/norm(w);
zs = (v'*w)*w;
end

function [X, Xs, Y, Ys] = child(t);
%Child
X  =  5*cos(t);  Y  =  5*sin(t);
Xs = -5*sin(t);  Ys =  5*cos(t);
end

[Andrea Scodro]

Obj3 is connected to Obj2 not to Obj1.

For the Obj3 I need to modify the f function. The new path is given by the position of the Obj2 along time.

function zs = f(t,z)
%
[X Xs Y Ys] = child(t);% <---For the second interaction I need to edit this line with the path of Obj2
v =[Xs; Ys];
w =[X-z(1); Y-z(2)];
w = w/norm(w);
zs = (v'*w)*w;
end

Another way to see the problem is that in the solution the distance between Obj2 and Obj3 is always 3 during the kinematics because they are connected by a rigid bar. If I do as you told me at the end of the kinematics the positions of Objs is:

Obj1:( 0 , 5 )

Obj2: ( 2.5184 , 3.3637 )

Obj3: ( 3.6499 , 0.2347 )

The distance between Obj1 and Obj2 is 3 and it is correct, but the distance between Obj2 and Obj3 result 3.33 more or less and it is uncorrect.

The solution might be right if Obj3 is connected to Obj1 but this is not the case.

I'm not familiar with MATLAB but I think I should interpolate (maybe with a spline) the result obtained for Obj2 and write a new function similar to child(t):

function [X, Xs, Y, Ys] = Obj2(t);
%Obj2
X  =  y(:,1);  Y  =  y(:,2); % Or something to interpolate y results
Xs = diff(X);  Ys =  diff(Y); % Something to calculate the differential of the X and Y functions over time
end

Now write a new function similar to f:

function zs = f2(t,z)
%
[X Xs Y Ys] = Obj2(t);
v =[Xs; Ys];
w =[X-z(1); Y-z(2)];
w = w/norm(w);
zs = (v'*w)*w;
end

And calculate the Obj3 path with the command:

% Object 3
y02 = [5; -6];
[t y2] = ode45(@f2,t,y02);

I think the idea is good but I don't know how to translate it into working code .

[Cris LaPierre]

The challenge is that ode solvers in MATLAB use a single time point at a time as they solve. For example, it solves for y(2,:) at time t(2) using t(1) and y0. It then solves for y(3,:) at t(3) using t(2) and y(2,:). At no point will your function f2 have access to all values of t and y2. Your solution will have to take that into account.

One way around this is to explicitly pass in the solution from obj1 as additional parameters to ode45. See this example for one way to do this.

Also, be sure to keep variable scope in mind. You use y inside your obj2 function without ever passing it in. If you need a refresher on functions, see this page.

[Cris LaPierre]

One more comment that, without an equation for obj3's position, you will likely discover approximation errors. You can certainly get closer to 3m for the distance between Obj2 and Obj3, but it will likely not be exactly 3.

[Cris LaPierre]

I think I may have found the same code as you online: The Tractrix and Similar Curves

I had a little more time, so here's a potential solution. This uses the numerical approach you proposed. As a caution, numerical methods are going to have approximation errors, so the distance of your linkage between Obj2 and Obj3 will not be exactly 3.

The second call to ode45 passes the result for Obj2 as additional inputs, making them avaiable to the odefunction f2.

y0 = [5; -3];
y02 = [5; -6];
% By specifying start and end points, I let ode45 choose time steps
tspan = [0 pi/2];
[t y] = ode45(@f1,tspan,y0); % Obj2
% Solve at same time points as Obj2 (makes it easier to calculate distance)
[t2 y2] = ode45(@f2,t,y02,[],t,y); % Obj3

[X, Xs, Y, Ys] = child1(t);
plot(X,Y,':','DisplayName','Obj1')
hold on;
plot(y(:,1),y(:,2),'DisplayName','Obj2');
plot(y2(:,1),y2(:,2),'DisplayName','Obj3');
hold off;
axis([-6 6 -6 6]);
axis('square');
legend

Calculate distances between objects

% Obj 1-2
max(sqrt(sum(([X,Y]-y(:,1:2)).^2,2)))

ans = 3.0000

% Obj 2-3
max(sqrt(sum((y(:,1:2)-y2).^2,2)))

ans = 3.0587

The solution can be improved by decreasing the step size.

tspan = [0:0.001:pi/2];
[t y] = ode45(@f1,tspan,y0); % Obj2
[t2 y2] = ode45(@f2,t,y02,[],t,y); % Obj3
max(sqrt(sum((y(:,1:2)-y2).^2,2)))

ans = 3.0020

Local Functions

function zs = f1(t,z)
[X Xs Y Ys] = child1(t);
v =[Xs; Ys];
w =[X-z(1); Y-z(2)];
w = w/norm(w);
zs = (v'*w)*w;
end
function [X, Xs, Y, Ys] = child1(t)
X  =  5*cos(t);  Y  =  5*sin(t);
Xs = -5*sin(t);  Ys =  5*cos(t);
end

Below are the modified functions to numerically solve for the velocity of Obj3, which ode45 then integrates to X,Y position data. Both functions have the additional inputs t1,y1, the solution from ode45 for Obj2.

function zs = f2(t,z,t1,y1)
% Inputs t1 and y1 are the ode45 solution for obj2
[X Xs Y Ys] = child2(t,t1,y1);
v =[Xs; Ys];
w =[X-z(1); Y-z(2)];
w = w/norm(w);
zs = (v'*w)*w;
end
function [X, Xs, Y, Ys] = child2(t,t1,y1)
% Inputs t1 and y1 are the ode45 solution for obj2
% Position is straight forward interpolation to solver time t
% Velocity is (change in position)/(change in time) interpolated to solver time t
X  =  interp1(t1,y1(:,1),t,"spline");  Y  =  interp1(t1,y1(:,2),t,"spline");
Xs =  interp1(t1,[0; diff(y1(:,1))./diff(t1)],t,"spline");  
Ys =  interp1(t1,[0; diff(y1(:,2))./diff(t1)],t,"spline");
end

[Andrea Scodro]

Amazing!! This is exactly what I was looking for!! Thanks a lot my friend.