Asked by Eliraz Nahum
on 28 Sep 2018

Hello everyone, I have a problem of 2 coupled second order equations of the form that is described in the photo I attached. By the way – W(t) and T(t) are quite simple functions, and I plan calling another external function when using it. The reason I don't include them above is the fact that they have a "steps" form and hence time-interval dependent. As I see the only option is a numeric solution by using the ode45 function. I want to make sure that using the above function is possible at all by doing the manipulation that is also described in the attached photo. I will be happy to hear what is your advice – is there another good way to solve this problem? If not, is the methodic way described above should work? How should I write the command when using the ode45 function? [t,q]=ode45(@(t,q) F(t,q),t_span,q0) While F is a 2x1 vector that consists each of the above eauations.

Answer by Bruno Luong
on 30 Sep 2018

Edited by Bruno Luong
on 30 Sep 2018

Accepted Answer

Bellow is the code with 2 methods.

Observation: piecewise RUNGE-KUTTA method is more accurate whenever the trajectory is orbital around the attractor (0,0)

% S(t) := T(t) / W(t),

% S(t) = s(i) for tbin(i) <= t < tbin(i+1);

%

n = 10; % number of interval where S is constant

tbin = [0 cumsum(0.1+rand(1,n))];

s = rand(1,n);

s(end+1) = s(end);

q0 = rand(4,1);

% Cauchy condition

tspan = tbin([1 end]);

[t,q] = ode15s(@(t,q) odefun(t,q,tbin,s), tspan, q0);

x = q(:,1);

y = q(:,3);

plotsol(1,t,x,y,tbin,'single run ode15s');

%

t = [];

q = [];

for k=1:n

subts = tbin([k k+1]);

[tk,qk] = ode45(@(t,q) odefun(t,q,tbin,s), subts, q0);

t = [t; tk];

q = [q; qk];

q0 = q(end,:);

end

x = q(:,1);

y = q(:,3);

plotsol(2,t,x,y,tbin,'piecewise ode45');

function plotsol(fignum,t,x,y,tbin,str)

% plot the solution

fig = figure(fignum);

clf(fig);

% make figure wider

p = get(fig,'position');

p(3) = 800;

set(fig,'position',p);

subplot(1,2,1),

h = plot(t,[x,y]);

xlabel('t');

hold on

n = length(tbin);

for k=2:n

plot(tbin(k)+[0 0],ylim(gca),'-.k');

end

legend(h,'x','y');

title(str);

subplot(1,2,2),

plot(x,y);

xlabel('x');

ylabel('y');

title(str);

end

function qp = odefun(t,q,tbin,s)

x = q(1);

xp = q(2);

y = q(3);

yp = q(4);

[~,i] = histc(t,tbin);

S = s(i);

a = sqrt(x.^2+y.^2).^3;

b = sqrt(xp.^2+yp.^2);

c = S./b;

xpp = -x./a + xp.*c;

ypp = -y./a + yp.*c;

qp = [xp; xpp; yp; ypp];

end

Bruno Luong
on 30 Sep 2018

Try this:

if 1

t1=32.2;

w0=6169;

w11=2695;

t2=138.2;

w12=2082;

w21=397;

w22=100;

T1=26950;

T2=3973;

GM = 1.4077 * 10^16;

Time=1000;

g = 32.17;

tbin = [0 t1 t2 Time];

sfun1 = @(t) (g*T1)./interp1([0 t1],[w0 w11],t);

sfun2 = @(t) (g*T2)./interp1([t1 t2],[w12 w21],t);

sfun3 = @(t) zeros(size(t));

s = {sfun1, sfun2, sfun3, sfun3};

n = 3;

tspan = [0 Time];

v0 = 200;

gama = deg2rad(80);

%Starting conditions

x0 = 0;

y0 = 2.0926 * 10^7;

xdot0 = v0 * cos(gama);

ydot0 = v0 * sin(gama);

q0 = [x0; xdot0; y0; ydot0];

else

% Test case

n = 10;

GM = 1;

tbin = [0 cumsum(0.1+rand(1,n))];

s = rand(1,n);

s(end+1) = s(end);

q0 = rand(4,1);

tspan = tbin([1 end]);

end

%

t = [];

q = [];

for k=1:n

subts = tbin([k k+1]);

[tk,qk] = ode45(@(t,q) odefun(t,q,tbin,s,GM), subts, q0);

t = [t; tk];

q = [q; qk];

q0 = q(end,:);

end

x = q(:,1);

y = q(:,3);

plotsol(2,t,x,y,tbin,'piecewise ode45');

function plotsol(fignum,t,x,y,tbin,str)

% plot the solution

fig = figure(fignum);

clf(fig);

% make figure wider

p = get(fig,'position');

p(3) = 800;

set(fig,'position',p);

subplot(1,2,1),

h = plot(t,[x,y]);

xlabel('t');

hold on

n = length(tbin);

for k=2:n

plot(tbin(k)+[0 0],ylim(gca),'-.k');

end

legend(h,'x','y');

title(str);

subplot(1,2,2),

plot(x,y);

xlabel('x');

ylabel('y');

title(str);

end

function qp = odefun(t,q,tbin,s,GM)

x = q(1);

xp = q(2);

y = q(3);

yp = q(4);

[~,i] = histc(t,tbin);

if isnumeric(s)

S = s(i);

else

% function handle

S = feval(s{i},t);

end

a = sqrt(x.^2+y.^2).^3;

b = sqrt(xp.^2+yp.^2);

c = S./b;

xpp = -GM*x./a + xp.*c;

ypp = -GM*y./a + yp.*c;

qp = [xp; xpp; yp; ypp];

end

Jan
on 25 Oct 2018

@Bruno: [~,i] = histc(t,tbin); S=s(i) can cause discontinuities in the function to be integrated. Matlab's ODE integrators cannot handle this reliably, see http://www.mathworks.com/matlabcentral/answers/59582#answer_72047. ODE45 might find a trajectory, stop due to a too small stepsize or in the worst case calculate a "result" which is far beyond the actual bounds defined by the error tolerances.

In your case you limit the time span to one specific tbin([k k+1]), such that histc find the same bin during a run of the integration. So wouldn't it be cheaper to determine the contstant outside the integration instead of calling histc for each call of the function?

Bruno Luong
on 25 Oct 2018

Depend on type of discontinuity, as long as it does not "stiffness", affects only low-order terms, the error estimate and adaptive time step work just well, as in the case of this rocket eqt we discuss here.

I showed in my code it works just fine numerically by splitting intervals or integrate straight through in a single shot. ODE45 handles it just fine, and I don't even need to use ODE15s.

As for use TSPAN, yes, it can sync the break points, but I can't see how one can remove HISTC (which is not expansive either) since the ode function never knows before hand which interval ODE45 is working on. However one can loop in the intervals as I showed in my code for-loop that call sequentially ODE45, that's probably the cleanest way to handle discontinuity.

I and Walter went over this discussion already, see the comments part of the question.

Sign in to comment.

Answer by Eliraz Nahum
on 29 Sep 2018

Sorry guys but I still can't understand if the code I wrote os feasible...

Answer by Eliraz Nahum
on 1 Oct 2018

Thank you all for your answers. I read each one and learned a lot.

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 6 Comments

## Walter Roberson (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/421271-2-coupled-second-order-equations#comment_615796

## Bruno Luong (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/421271-2-coupled-second-order-equations#comment_615801

## Walter Roberson (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/421271-2-coupled-second-order-equations#comment_615900

## Bruno Luong (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/421271-2-coupled-second-order-equations#comment_615905

## Walter Roberson (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/421271-2-coupled-second-order-equations#comment_615919

## Bruno Luong (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/421271-2-coupled-second-order-equations#comment_615924

Sign in to comment.