# Solving of three non-linear simultaneous differential equations with given boundary conditions

I am trying to solve three simultaneous differential equations in Mathematica-10 but I am not able to get the solutions. After running the following command as I am pressing shift+enter it just simply shows the equation again:

My equations are:

{x1'[t]-w1[t]*l2*x2[t]-l1*(v-x1[t])==0,
x2'[t]+w1[t]*l1*x1[t]-l2*(v-x1[t])==0,
w1'[t]-k*x2[t]*(v-x1[t])==0}


with boundary conditions:

{x1[0]==0, x2[0]==0,w1[0]==15,x1[infinity]=v}


After solving these equations i need the values of x1, x2, and w1.

For solving these equations i ma running the following command:

DSolve[{
x1'[t]-w1[t]*l2*x2[t]-l1*(v-x1[t])==0,
x2'[t]+w1[t]*l1*x1[t]-l2*(v-x1[t])==0,
w1'[t]-k*x2[t]*(v-x1[t])==0,
x1[0]==0,x2[0]==0, w1[0]==50, x1[Infinity]==v
},{x1, x2, w1}, t]


Can anyone suggest me the correct answer of the above equations?

• Welcome to Mathematica.SE! I formatted your question. Please learn how to do it your self, and read this to see how to ask good questions. – yohbs Sep 24 '15 at 11:07
• DSolve will return an analytic solution. Are you sure such a solution exists? – march Sep 24 '15 at 17:06
• Unless this is meant to be an eigenvalue problem, there is one too many boundary conditions for the third-order system of ODEs. – bbgodfrey Sep 28 '15 at 2:12
• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. – bbgodfrey Sep 28 '15 at 14:07

As march suggested in a comment, it appears that DSolve cannot solve this problem. However, quite a good analytical approximation can be obtained, as presented below. First, however, it is useful to examine a numerical solution in order to gain insight. Parameters {l1 -> 1, l2 -> 2, k -> 1, v -> 1} are used here, although similar results are obtained for other choices, provided that the parameters are positive and somewhat less than w1[0].

s = NDSolve[{x1'[t] - w1[t]*l2*x2[t] - l1*(v - x1[t]) == 0,
x2'[t] + w1[t]*l1*x1[t] - l2*(v - x1[t]) == 0,
w1'[t] - k*x2[t]*(v - x1[t]) == 0, x1[0] == 0, x2[0] == 0,
w1[0] == 50} /. {l1 -> 1, l2 -> 2, k -> 1, v -> 1},
{x1[t], x2[t], w1[t]}, {t, 0, 5000}];


Multiple plots are needed to display the diverse scales present in the solution.

Plot[Evaluate[{x1[t], x2[t], w1[t]} /. s], {t, 0, 5000},
PlotRange -> All, AxesLabel -> {t, "x1, x2, w1"}]


Plot[Evaluate[{x1[t], x2[t], w1[t]} /. s], {t, 0, 5000},
PlotRange -> {Automatic, {-.4, 1.2}}, AxesLabel -> {t, "x1, x2, w1"}]


Plot[Evaluate[{x1[t], x2[t], w1[t]} /. s], {t, 0, 5},
PlotRange -> {Automatic, {-.1, .2}}, AxesLabel -> {t, "x1, x2, w1"}]


Two interesting features are immediately apparent. First, x1 is equal to v at large t, as desired by the OP, even without an imposed boundary condition there. (This occurs only for positive l1 and l2.) Second, the characteristic timescale for variation of x1 and x2 is far shorter than that for w1 before all three become essentially constant at late time. This suggests that reasonably accurate expressions for x1 and x2 can be obtained ignoring the variation in w1. For convenience, temporarily designate this quasi-constant quantity as w0 and solve for the remaining variables.

s0 = FullSimplify[ExpToTrig[Flatten@DSolve[
{x1'[t] - w0*l2*x2[t] - l1*(v - x1[t]) == 0,
x2'[t] + w0*l1*x1[t] - l2*(v - x1[t]) == 0,
x1[0] == 0, x2[0] == 0}, {x1[t], x2[t]}, t] /.
l1^2 - 4 l2^2 w0 - 4 l1 l2 w0^2 -> -c]]
(* {x1[t] -> (4 E^(-((l1 t)/2))
l2 v w0 (Sqrt[c] E^((l1 t)/2) l2 - Sqrt[c] l2 Cos[(Sqrt[c] t)/2] +
l1 (l2 + 2 l1 w0) Sin[(Sqrt[c] t)/2]))/(Sqrt[c] (c + l1^2)),
x2[t] -> (4 E^(-((l1 t)/2))
v w0 (-Sqrt[c] E^((l1 t)/2) l1^2 + Sqrt[c] l1^2 Cos[(Sqrt[c] t)/2] +
(l1^3 + 2 l2^3 + 2 l1 l2^2 w0) Sin[(Sqrt[c] t)/2]))/(Sqrt[c] (c + l1^2))} *)


c also is introduced for ease of visualization and manipulation. For w0 somewhat larger than the other parameters, c is approximately equal to 4 l1 l2 w0^2. Now, obtain limiting expressions for x1 and x2 for times long compared their transient behavior.

MapAt[Limit[#, t -> Infinity, Assumptions -> l1 > 0 && c > 0] &, s0, {All, 2}];
suv = Simplify[% /. c -> -(l1^2 - 4 l2^2 w0 - 4 l1 l2 w0^2)] /. w0 -> w1[t]
(* {x1[t] -> (l2 v)/(l2 + l1 w1[t]), x2[t] -> -((l1^2 v)/(l2^2 + l1 l2 w1[t]))} *)


We now can obtain w1, again using DSolve.

First@DSolve[Simplify[w1'[t] - k*x2[t]*(v - x1[t]) == 0 /. suv], w1[t], t]
(* {w1[t] -> InverseFunction[l2^2 Log[#1] + 2 l1 l2 #1 + (l1^2 #1^2)/2 &]
[-((k l1^3 t v^2)/l2) + C[1]]} *)


This would be a completely symbolic solution, except that C[1] can be obtained from w1[0] == 50 only after the other parameters have been given numerical values. For instance,

f = w1[t] /. % /. C[1] -> t0 /. {l1 -> 1, l2 -> 2, k -> 1, v -> 1};
FindRoot[(f /. t -> 0) == 50, {t0, 1500}]
(* {t0 -> 1465.65} *)


with C[1] replaced by t0 to keep FindRoot happy.

Alternatively, and more simply,

sw = First@DSolve[{Simplify[w1'[t] - k*x2[t]*(v - x1[t]) == 0 /. suv /.
{l1 -> 1, l2 -> 2, k -> 1, v -> 1}], w1[0] == 50}, w1[t], t]
(* {w1[t] -> InverseFunction[4 Log[#1] + 4 #1 + #1^2/2 &]
[-(t/2) + 2 (725 + 2 Log[50])]} *)


Altogether, the mostly symbolic solution for this set of parameters is

{x1[t], x2[t], w1[t]} /. suv /. sw /. {l1 -> 1, l2 -> 2, k -> 1, v -> 1}


and plotting it, for instance by

Plot[Evaluate[{x1[t], x2[t], w1[t]} /. suv /. sw /. {l1 -> 1, l2 -> 2, k -> 1, v -> 1}],
{t, 0, 5000}, PlotRange -> {{0, 5000}, {-.4, 1.2}}, WorkingPrecision -> 30]


yields a curve indistinguishable to the eye from the second curve above.

The results above are entirely symbolic except for the determination of C[1], which employes FindRoot. This constant can be obtained symbolically as well. From the earlier analysis,

w1[t] /. First@DSolve[Simplify[w1'[t] - k*x2[t]*(v - x1[t]) == 0 /. suv], w1[t], t]


gives w1[t] as a function of C[1] - ((k l1^3 t v^2)/l2). Consequently, C[1] - ((k l1^3 t v^2)/l2) can be expressed as a function of w1[t] by

InverseFunction[Head[%]]
(* l2^2 Log[#1] + 2 l1 l2 #1 + (l1^2 #1^2)/2 & *)


which at w1[0] = 50 yields C[1]

%[50]
(* 1250 l1^2 + 100 l1 l2 + l2^2 Log[50] *)


For the set of parameters, {l1 -> 1, l2 -> 2, k -> 1, v -> 1}, this equals 1450 + 4 Log[50], or 1465.65, as computed previously.