Solving a system of nonlinear coupled ODEs

Question

I tried using NDSolve to solve the following system of equations:

t x'[t] == -x[t] + y[t], 
t y'[t] == -5 t^2/x[t]^2 + x[t] - y[t],       
x[1] == 4, x[100] == 1

It's weird. The system tells me there is infinity at the boundary $t=1$, so I change the boundary from $t=1$ to $t=2$, and I get the same message again.

If I eliminate the nonlinear term -5 t^2/x[t]^2, the function can be solved analytically. So I do not know whether these equations are well posed or not with such boundary conditions just for x[t]

Anyone have a suggestion?

Answer 1

As noted above, you need an initial condition for y. You can then solve the odes given x[1] == 4 and y[1] == y1 where y1 is some real number. By solving the equations for different values of y1, you can see if there is a value such that x[100] == 1. In Version 9, it turns out the ParametricNDSolve will help you do this. Make y1 a parameter:

{xsoln, ysoln} = {x, y} /. ParametricNDSolve[{
  t x'[t] == -x[t] + y[t], 
  t y'[t] == -5 t^2/x[t]^2 + x[t] - y[t],
  x[1] == 4,
  y[1] == y1
}, {x, y}, {t, 1, 100}, {y1}]

Now locate a region where y1 is likely to give the desired result:

Plot[xsoln[y1][100], {y1, 43.35, 43.5}, GridLines -> {{}, {1}}]

Next use FindRoot to get a precise value and call it y1star:

y1star = y1 /. FindRoot[xsoln[y1][100] == 1, {y1, 43.4}]

Finally, confirm that it produces the desired result

Plot[{xsoln[y1star][t], ysoln[y1star][t]}, {t, 1, 100}, 
  PlotRange -> {-10, All}, GridLines -> {{}, {1}}]

Answer 2

Although the solution by Nasser is entirely satisfactory, this alternative approach may be of interest. Designate the two equations in the question as eq1 and eq2, respectively. Then, since y has no boundary conditions, eliminate it from the system.

Flatten@Solve[Subtract @@ eq1 + Subtract @@ eq2 == 0, y'[t]]
eq3 = Simplify[D[eq1, t] /. %]
(* {y'[t] -> (-5 t - x[t]^2 x'[t])/x[t]^2} *)
(* (5 t)/x[t]^2 + 3 x'[t] + t x''[t] == 0 *)

Then, solve the resulting equation for x, and from that result determine y as well (from eq1).

{sx, sxp} = NDSolveValue[{eq3, x[1] == 4, x[100] == 1}, {x, x'}, {t, 1, 100}, 
    Method -> {"Shooting", "StartingInitialConditions" -> x'[1] == 39.5}];
Plot[{sx[t], sx[t] + t sxp[t]}, {t, 1, 100}, PlotRange -> {-20, All}]

Answer 3

As pointed out by Nasser, you need to have an initial/boundary condition for y(t). Here I assumed a random condition on y(t)

sol = NDSolve[{t x'[t] == -x[t] + y[t], t y'[t] == -5 t^2/x[t]^2 + x[t] - y[t], 
x[1] == 4, y[1] == 1}, {x, y}, {t, 1, 2}]

Plot[{x[t] /. sol, y[t] /. sol}, {t, 1, 2}, PlotRange -> All, 
PlotStyle -> {Blue, Red}, PlotLegends -> Placed[{"x[t]", "y[t]"}, Above]]