# How to solve two-point boundary nonlinear ordinary differential equations （ODEs）with free final time numerically?

I am trying to solve two-point boundary nonlinear ordinary differential equations (ODEs) with free final time. I have failed many times to apply Mathematica to solve this problem and would really appreciate if someone can give me an answer.

Given variables: v=1; w=0.00003; r=0.0002; {xo,yo} = {20,20}:

The ODEs are:

x'[t] ==vCos[q[t]]+wx[t]+ry[t]+(x[t]-xo)/((x[t]-xo)^2+(y[t]-yo)^2)^(3/2);
y'[t] ==vSin[q[t]]-rx[t]+wy[t]+(y[t]-yo)/((x[t]-xo)^2+(y[t]-yo)^2)^(3/2);
q'[t] ==-r+3(Cos[q[t]^2]-Sin[q[t]^2])(x[t]-xo)(y[t]-yo)/((x[t]-xo)^2+(y[t]-yo)^2)^(5/2)+
3Cos[q[t]]Sin[q[t]]((y[t]-yo)^2-(x[t]-xo)^2)/((x[t]-xo)^2+(y[t]-yo)^2)^(5/2);


where x'[t], y'[t] are the known motion dynamics of the system and q'[t] is the evolution of control law which is already known too. Given an initial known condition [x[t=0], y[t=0]]=[0,0] and the final target condition {x[t=tf], y[t=tf]} = {100,100} where the final time tf is unknown, the question is how to find the initial control input q[t=0] which makes the system situation changes from the initial condition {0,0} to the final target condition {100,100} with the minimal time tf?

• Please provide the actual equations, so that readers can experiment with them. Be sure to write them in Mathematica format. Sep 21, 2016 at 0:10
• Thanks for your reply and I have given the detailed equations, please help me on solving the problem. Looking forwards to hearing from you. Sep 21, 2016 at 8:54
• Is vCos a free parameter or should it be v*Cos[q[t]]? Note Cos(x) means Cos times x, which is almost never what is wanted. Sep 21, 2016 at 9:47
• You are right that it represents v*Cos[q[t]], please help me on solving this problem. Sep 21, 2016 at 10:17
• I have tried the shooting method to solve this problem where the initial control input 'q[t=0] ' ranges from [0 2*pi]. However, I don't know how to formulate the problem to Mathematica form when applying the shooting method as the final time 'tf' is unknown. @ bbgodfrey @ Michael E2 Sep 21, 2016 at 14:17

Set up the integration as follows:

eq = {x'[t] == v Cos[q[t]] + w x[t] +
r y[t] + (x[t] - xo)/((x[t] - xo)^2 + (y[t] - yo)^2)^(3/2),
y'[t] == v Sin[q[t]] - r x[t] +
w y[t] + (y[t] - yo)/((x[t] - xo)^2 + (y[t] - yo)^2)^(3/2),
q'[t] == -r + 3 (Cos[q[t]^2] - Sin[q[t]^2]) (x[t] - xo) (y[t] - yo)/((x[t] - xo)^2 +
(y[t] - yo)^2)^(5/2) + 3 Cos[q[t]] Sin[q[t]] ((y[t] - yo)^2 -
(x[t] - xo)^2)/((x[t] - xo)^2 + (y[t] - yo)^2)^(5/2)}

tmax = 300;
s = ParametricNDSolveValue[{eq /. {v -> 1, w -> 3 10^-5, r -> 2 10^-4, xo -> 20, yo -> 20},
x[0] == 0, y[0] == 0, q[0] == q0}, {x, y, q}, {t, 0, tmax}, {q0}];


sol = FindRoot[(Through[s[q0][tf]] - 100)[[1 ;; 2]], {tf, 140}, {q0, .88}, Evaluated -> False]
(* {tf -> 141.238, q0 -> 0.885835} *)


I obtained the initial guesses by plotting s[q0] for various values of test, for instance.

test = 140;
Plot[{Through[s[q0][test]][[1]] - 100, Through[s[q0][test]][[2]] - 100}, {q0, 0.8, 1}]


For completeness, here is a Plot of the integrated functions with parameters that satisfy the desired boundary conditions.

Plot[Evaluate@Through[s[q0][t] /. sol], {t, 0, tf /. sol}, AxesLabel -> {t, "x, y, q"}]


In response to comments below, a slightly more efficient approach is

s = ParametricNDSolveValue[{eq /. {v -> 1, w -> 3 10^-5, r -> 2 10^-4, xo -> 20, yo -> 20},
x[0] == 0, y[0] == 0, q[0] == q0}, {x, y, q}, {t, 0, tf}, {q0, tf}];
sol = Quiet@FindRoot[(Through[s[q0, tf][tf]] - 100)[[1 ;; 2]],
{tf, 140}, {q0, 1.0}, Evaluated -> False]
(* {tf -> 141.238, q0 -> 0.885835} *)


as before. (Quiet has been added to suppress unimportant error messages produced when guesses are not so good.) Interestingly, a second solution also exists.

sol = Quiet@FindRoot[(Through[s[q0, tf][tf]] - 100)[[1 ;; 2]],
{tf, 140}, {q0, 0.6}, Evaluated -> False]
(* {tf -> 141.379, q0 -> 0.711415} *)


tf is slightly larger for this second solution.

• Thanks very much for your contribution. What does Evaluated -> False mean and the purpose of the usage of tmax? @bbgodfrey Sep 21, 2016 at 14:51
•  Evaluated -> False prevents FindRoot from trying to evaluate s without inserting numerical values. tmax is just an upper bound on the integration that I chose to be large enough that tf was likely to be less than tmax. Sep 21, 2016 at 14:57
• Great appreciation for your help. In that case, I need to choose a large engough tmax when I solve the ODEs with different initial and final conditions. Thanks again for your code and I will have a try on it for different initial and final conditions. @ bbgodfrey Sep 21, 2016 at 15:03
• @JohnBai In response to your last comment, I modified the code slightly to eliminate the need to guess tmax. In the process, I noticed a second solution at smaller q0. It has a slightly larger value of tf, so the original solution still is the one you wanted. However, you need to be aware that there may be two solutions when running different parameters. Sep 21, 2016 at 19:41
• Thanks very much for your help. I tried to use the code on different initial and final conditions. It works well except that for each case the initial guess on test and q0 should be conducted. Thanks for your reply. @bbgodfrey Sep 22, 2016 at 7:58