2 coupled second order equations

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.

[Walter Roberson]

If you have the Symbolic Toolbox then I recommend using the flow documented in the first example of odeFunction() to set up the equations.

However, step functions cannot be used with any of the ode*() routines, all of which assume continuous derivatives. Perhaps you should be using delay differential functions.

[Bruno Luong]

Uh uh, he told about step in the coefficients T/W, it has nothing to do with the smoothness of the solution, at least not directly at the same differential order.

If the non linearity is "nice", the solution might be stiff but solvable with some stiff method.

Secondly, nothing prevent user to solve the ode (using smooth solver) successively on each time interval on which the coefficients have no jump.

[Walter Roberson]

If the non linearity is "nice", the solution might be stiff but solvable with some stiff method.

Not if it crosses any step boundaries. If it does then it is not solvable using the ode*() routines, except perhaps as numeric nonsense.

Secondly, nothing prevent user to solve the ode (using smooth solver) successively on each time interval on which the coefficients have no jump.

Correct, but not something I felt like getting into at 6:30 AM after having been awake all night.

[Bruno Luong]

Not if it crosses any step boundaries. If it does then it is not solvable using the ode*() routines, except perhaps as numeric nonsense.

Wrong, it makes perfect sense (ever heard of "weak solution" or are you still thinking about classical derivative?), the ODE solver adapt the step, and detect the stiffness. I can cross the boundary with the small step, and the error is still small and ensure the stability.

Try this simple example to revise your opinion:

ode45(@(t,y) double(t>=1), [0 2], 0)

here I don't even need stiff solver for step RHS.

[Walter Roberson]

Runge–Kutta methods to estimate slopes are only valid if the slopes have continuous derivatives, which is not the case for step functions.

ode45() does not always detect the discontinuity, which just means that the answers produced are not always valid.

[Bruno Luong]

The non-stiff method such as Runge-Kutta is less accurate at the breakpoints, but it is a valid solution.

I don't know if you see the graph of the example I gives, but yes it detect the discontinuity indirectly by estimate the error term and it reduces the step when it cross the break point. The solution looks OK (zoom around t=1, you'll see it negotiates the break point very very well). The theoretical solution at final time t=2 is y=1 theorically, and ode45 has an error of 3.148e-5. Tell me if that accuracy does not satisfy you.

The ODE solver are more sophisticated than you might think.

Nothing telling me (beside you of course) that the solution is not valid.

Answer 1

[Bruno Luong]

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

[Eliraz Nahum]

hey everyone and thanks for your kind help. I am sitting helpless, but thanks to you I am better now

while waiting to your (good) and sophisticated answer, my friend and I tried to simplify the problem. We got 2 different answers (the x values in the graph) and can't figure out where is the problem. we have to hand it out today, and really don't know what to do. I attached the original physical problem that we were presented and our 2 solutions. we will be glad to know what is true and why the other one is wrong.

1st code:

clear all
close all
clc

t1=32.2;
w0=6169;
w11=2695;
t2=138.2;
w12=2082;
w21=397;
w22=100;
T1=26950;
T2=3973;
x0=0;
y0=2.0926*(10^7);
xdot0=34.7;
ydot0=196.9;

g=32.15;
GM=1.4077*(10^16); 

t_delta=1;

counter=1;
time=0;
v_time=[];

x=[];
y=[];
xdot=[];
ydot=[];
x2dot=[];
y2dot=[];

x(counter)=x0;
y(counter)=y0;
xdot(counter)=xdot0;
ydot(counter)=ydot0;
division1=T1/w0;
x2dot(counter)=(-GM*(x(1)/((x(1)^2 + y(1)^2)^1.5)))+(g*division1)*(xdot(1)/((xdot(1)^2 + ydot(1)^2)^0.5));
y2dot(counter)=(-GM*(y(1)/((x(1)^2 + y(1)^2)^1.5)))+(g*division1)*(ydot(1)/((xdot(1)^2 + ydot(1)^2)^0.5));

r=sqrt(x0^2+y0^2);
size=r;

while (size>=r)
time=time+t_delta;

if  time<t1
    W=w0+(((w11-w0)/t1)*time);
    T=T1;
elseif (time>=t1)&&(time<t2)
    W=w12+(((w21-w12)/(t2-t1))*(time-t1));
    T=T2;
else
    W=w22;
    T=0;
end
TvsW=T/W;

counter=counter+1;

x(counter)=x(counter-1)+xdot(counter-1)*t_delta;
y(counter)=y(counter-1)+ydot(counter-1)*t_delta;

xdot(counter)=xdot(counter-1)+x2dot(counter-1)*t_delta;
ydot(counter)=ydot(counter-1)+y2dot(counter-1)*t_delta;

x2dot(counter)=(-GM*(x(counter)/((x(counter)^2 + y(counter)^2)^1.5)))+(g*TvsW)*(xdot(counter)/((xdot(counter)^2 + ydot(counter)^2)^0.5));
y2dot(counter)=(-GM*(y(counter)/((x(counter)^2 + y(counter)^2)^1.5)))+(g*TvsW)*(ydot(counter)/((xdot(counter)^2 + ydot(counter)^2)^0.5));

v_time(counter)=time;

size=sqrt(x(counter)^2+y(counter)^2);
end

plot(x,y,'.b')
xlim([0 (10*10^6)])
ylim([0 (3*10^7)])
figure
plot(v_time,x,'.r')
hold on
plot(v_time,y,'.g')
figure
plot(v_time,xdot,'.r')
hold on
plot(v_time,ydot,'.g')

figure(1)
hold on
for counter=0:0.001*pi:2*pi
    plot(r*cos(counter),r*sin(counter),'.m')
    hold on
end
hold off

the second code:

close all;
clear all;
clc;

%Basic data
GM = 1.4077 * 10^16;
v0 = 200;
gama = deg2rad(80);
g = 32.17;
dt=1;
Time=100;

%Starting conditions
x0 = 0;
y0 = 2.0926 * 10^7;
xdot0 = v0 * cos(gama);
ydot0 = v0 * sin(gama);

for n=1:Time

      T = 26950;
      w = 6169;
      t=0;

      s1(1,n) = 30.59 + rand * 2 * 1.61;
      s2(1,n) = 131.29 + rand * 2* 6.91;

      x_val(1,1) = x0;
      x_dot_val(1,1) = xdot0; %=x1dot

      y_val(1,1) = y0; % =R world
      y_dot_val(1,1) = ydot0; %=y1dot

      while sqrt(x_val ^ 2 + y_val ^ 2) >= y0   

          v = ( x_dot_val ^ 2 + y_dot_val ^ 2 ) ^0.5;
          aT = ( g * T ) / w;

          x_2dot_val = ( -GM * x_val) / ((x_val ^ 2 + y_val ^ 2)^1.5) + (aT * x_dot_val) / v;
          y_2dot_val = ( -GM * y_val) / ((x_val ^ 2 + y_val ^ 2)^1.5) + (aT * y_dot_val) / v;

          x_val = x_val + dt * x_dot_val;
          x_dot_val = x_dot_val + dt * x_2dot_val;

          y_val = y_val + dt * y_dot_val;  
          y_dot_val = y_dot_val + dt * y_2dot_val;   

              if t<=s1(1,n)
                  w = ((2695-6169)/(32.2-0))*dt + 6169;
              elseif t>s2(1,n)
                  w = 100;
                  T = 0;
                else
                  w = ((397 - 2082) / (138.2 - 32.2)) * dt + 2082;
                  T = 3973;
              end

      t=t+dt;

      x_vector(1,t) = x_val;
      y_vector(1,t) = y_val;
      x_dot_vector(1,t) = x_val;
      y_dot_vector(1,t) = y_val;
      x_2dot_vector(1,t) = x_val;
      y_2dot_vector(1,t) = y_val;
      end

      length_arc(1,n) = atan (x_val/y_val) * y0; %חישוב הקשת במעגל
      figure(1);
      plot(x_vector,y_vector);
      hold on;

end

% g = viscircles(x,y,0,0,y0);
% plot(g);
th = 0:pi/50:2*pi;
xunit = y0 * cos(th) ;
yunit = y0 * sin(th);
h = plot(xunit, yunit);
hold off
figure(2);
hist(length_arc,20);

[Bruno Luong]

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]

@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]

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.

Answer 2

[Eliraz Nahum]

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

[Bruno Luong]

It seems perfectly OK to me.

Answer 3

[Eliraz Nahum]

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