# Non-linear system of ODEs with positive functions

I would like to solve this non-linear system of ODEs, but I need to tell Mathematica to consider y[t] and x[t] to be always positive. How can I do?

   system = {
t*y'[t] + 3*y[t] == -y[t]/Sqrt[x[t]+y[t]],
t*x'[t] + 4*x[t] == y[t]/Sqrt[x[t]+y[t]],
y[1]==1,
x[1]==0
};
s = NDSolve[system, {x, y}, {a, 1, 1000}]


Update

I've tried the method of Alexei Boulbitch and it works fine.

However, if I slightly modify the equations above by multiplying both the r.h.s with a constant const=10^9, I obtain the following error message:

I don't understand what's the problem.

The same issue comes out if I choose a very high initial condition for y[1], such as y[1]=10^32 or higher.

I would like to include both these add-ons at once, how can I handle these problems?

• What do you mean by "tell Mathematica to consider y[t] and x[t] to be always positive"? If you want NDSolve to stop when x and y is almost 0, check document of WhenEvent. Aug 22 '19 at 10:30

Try the following: First, let us define your system as you did, but shift the initial condition a bit from zero. Otherwise, it gives rise to a singularity:

system = {t*y'[t] + 3*y[t] == -y[t]/Sqrt[x[t] + y[t]],
t*x'[t] + 4*x[t] == y[t]/Sqrt[x[t] + y[t]], y[1] == 1,
x[1] == 0.001};


Now, let us replace variables such that x(t)=φ[t]^2 and y[t]=ψ[t]^2. Like this, the variables x and yare always non-negative. And then let us rewrite the initial conditions in terms of the new variables, φ[t] and ψ[t]:

system2 = ReplacePart[system /. {x -> (φ[#]^2 &),
y -> (ψ[#]^2 &)}, {{3, 1} -> ψ[1], {4} -> φ[1] == Sqrt[0.001]}]


Now, let us solve it:

{nds1, nds2} = NDSolveValue[system2, {φ, ψ}, {t, 1, 10}]


and plot x(t)=φ[t]^2 and y[t]=ψ[t]^2:

 Plot[{nds1[t]^2, nds2[t]^2}, {t, 1, 10}, PlotStyle -> {Red, Blue},
PlotRange -> {0, 0.12}, AxesLabel -> {Style["t", Italic, 16],
Row[{Style["x(t), ", Red, Italic, 16], Style["y(t)", Blue, Italic, 16]}]}]


which returns the following:

Another possibility would be to replace x and y by exponents, as follows:

system2 =
ReplacePart[system /. {x -> (Exp[φ[#]] &),
y -> (Exp[ψ[#]] &)}, {{3} -> ψ[1] == 0, {4} -> φ[1] == Log[0.001]}]


The rest is exactly as in the previous case, and the resulting plot

Plot[{Exp[nds1[t]], Exp[nds2[t]]}, {t, 1, 10},
PlotStyle -> {Red, Blue}, PlotRange -> {0, 0.12},
AxesLabel -> {Style["t", Italic, 16],
Row[{Style["x(t), ", Red, Italic, 16],
Style["y(t)", Blue, Italic, 16]}]}]


looks the same.

Have fun!

• Thank you so much! Why in this line, system2 = ReplacePart[system /. {x -> (φ[#]^2 &), y -> (ψ[#]^2 &)}, {{3, 1} -> ψ[1], {4} -> φ[1] == Sqrt[0.001]}], do you only set φ[1] and why do you use Sqrt[]? Shouldn't it be φ[1] == 0.001^2 ? Indeed, I have tried to put a very high initial condition for x[1]==10^16, but it gives the error "NDSolveValue::ndsz: At a == 1., step size is effectively zero; singularity or stiff system suspected." What could be the reason?
– Lele
Aug 22 '19 at 14:07
• Shouldn't it be φ[1] == 0.001^2 ? No, since here x=φ^2. Therefore, φ[1]=Sqrt[x[1]]. But in fact, you are free to invent any nonlinear transformation of your liking, such that it keeps x and y positive. Aug 22 '19 at 15:20
• Oh yes, sorry... But what about my second question? It doesn't seem to like very high initial conditions
– Lele
Aug 22 '19 at 15:30
• It is most probably that the initial condition x[1]==10^16 is much too high. I am not sure, but it seems that this value is on the limit admitted by Mma. I do not see, why do you need that high initial value. Are you sure that it is for x, rather than for y? If you believe that one of these functions contains a singularity, a good idea might be to separate the singularity out explicitly. Aug 23 '19 at 7:13
• Sorry, you're right, I meant y[1], but the problem is there anyway. For simplicity, I have posted a dimensionless version of the ODE I am interested in. I have tried to generalize it by including the huge physical constants. For example, the two add-ons are: 1) A very high initial condition for y[1]; 2) a multiplicative constant const=2*10^9 in front of the two y[t]/Sqrt[x[t] + y[t]]. If I include only the second condition, after NDSolveValue I get NDSolveValue::ndsz: At a == 1., step size is effectively zero; singularity or stiff system suspected`
– Lele
Aug 23 '19 at 7:39