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.
Addendum
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.
DSolve
will return an analytic solution. Are you sure such a solution exists? $\endgroup$