Solve Differential Equations with Boundary Condition at Variable Point

Question

I have a set of differential equations:

 w1'[u] == -1 - 1/w1[u] + 1/w2[u] +  1/w3[u],
 w2'[u] == -1 -  1/w1[u] +1/w2[u] +    1/w3[u],
 w3'[u] == -1 - 1/w1[u] + 1/w2[u] +  1/w3[u],
 w2R'[u] == -1 - 1/w2R[u] + 1/w2L[u] +  1/w4C[u],
 w2L'[u] == -1 -  1/w2R[u] +1/w2L[u] +    1/w4C[u],
 w4C'[u] == -1 - 1/w2R[u] + 1/w2L[u] +  1/w4C[u],

and in addition to some some normal boundary conditions

w1[Log[91.1876]] == 1/0.016887, w2[Log[91.1876]] == 1/0.03322, 
w3[Log[91.1876]] == 1/0.12

the boundary conditions for w2R,w2L,w4C are a bit more variable. Firstly, they should meet at a point $M$, which should be computed and is a priori unknown:

w2R[M]==w2L[M]==w4C[M]

and secondly at some other, a priori unknown point M2 the solutions should fulfil

 w1[M2]== w2R[M2] +w4C[M2]
 w2[M2]== w2R[M2]
 w3[M2]== w4C[M2]+3

How I can get the solutions of this system, i.e. M,M2, w1, w2,w3,w2R,w2L,w4C?

The system for w1,w2,w3 can be solved individually with NSolve, but my problem is how to solve the system for w2R,w2L,w4C, with the boundary conditions dependent on M and M2 that should be computed by Mathematica, too.

For a simpler set of differential equations I could simply solve, for example

w1'[u] == -1 ,
 w2'[u] == -1 ,
 w3'[u] == -1,

w1[Log[91.1876]] == 1/0.016887, w2[Log[91.1876]] == 1/0.03322, 
w3[Log[91.1876]] == 1/0.12

analytically using DSolve and likewise

 w2R'[u] == -1,
 w2L'[u] == -1 ,
 w4C'[u] == -1 ,

with DSolve without boundary conditions. Then the analytic solutions for w2R,w2L and w4C depend on three constants C1, C2, C3. Then I can compute the unknowns C1, C2, C3, MI and MU using ordinary NSolve and the boundary conditions

 w2R[M]==w2L[M]
 w2R[M]==w4C[M]
 w1[M2]== w2R[M2] +w4C[M2]
 w2[M2]== w2R[M2]
 w3[M2]== w4C[M2]+3

Unfortunately, this approach fails for more complicated differential equations that can't be solved analytically.

Answer 1

Because the right sides of the ODEs for {w2R, w2L, w4C} are equal, these three dependent variables are equal up to constants. Moreover, the boundary condition w2R[M]==w2L[M]==w4C[M] forces those constants to be zero. Thus, {w2R, w2L, w4C} are equal, independent of the value of M.

Similarly, because the right sides of the ODEs for {w1, w2, w3} are equal, these three dependent variables also are equal up to constants, which can be computed from boundary conditions.

c2 = (w2[u] - w1[u]) /. u -> Log[91.1876] 
    /. {w1[Log[91.1876]] -> 1/0.016887, w2[Log[91.1876]] -> 1/0.03322}
(* -29.1148 *)
c3 = (w3[u] - w1[u]) /. u -> Log[91.1876] 
    /. {w1[Log[91.1876]] -> 1/0.016887, w3[Log[91.1876]] -> 1/0.12}
(* -50.8838 *)

Hence, the six ODEs are reduced to two.

eqw = w1'[u] == (-1 - 1/w1[u] + 1/w2[u] + 1/w3[u]) 
    /. {w2[u] -> w1[u] + c2, w3[u] -> w1[u] + c3} 
(* w1'[u] == -1 + 1/(-50.8838 + w1[u]) + 1/(-29.1148 + w1[u]) - 1/w1[u] *)
eqLRC = w4C'[u] == (-1 - 1/w2R[u] + 1/w2L[u] + 1/w4C[u]) 
    /. {w2R[u] -> w4C[u], w2L[u] -> w4C[u]}
(*  w4C'[u] == -1 + 1/w4C[u] *)

Solutions to both can be obtained using DSolve. However, first consider the remaining boundary conditions, partially evaluated using the results obtained above.

{w1[M2] == w2R[M2] + w4C[M2], w2[M2] == w2R[M2], w3[M2] == w4C[M2] + 3} 
  /. {w2R[M2] -> w4C[M2], w2L[M2] -> w4C[M2], w2[M2] -> w1[M2] + c2, w3[M2] -> w1[M2] + c3}
(* {w1[M2] == 2 w4C[M2], -29.1148 + w1[M2] == w4C[M2], -50.8838 + w1[M2] == 3 + w4C[M2]} *)

With three equations and two unknowns, there is no solution. Until this inconsistency is resolved, little more can be done.

Addendum

For completeness, eqw can be solved as

Quiet@DSolve[{eqw, w1[Log[91.1876]] == 1/0.016887}, w1, u]

which returns only partly solved.

(* Solve[3.36112*10^-6 (0.0420198 Log[-51.909 + w1[u]] + 0.0342149 Log[-30.0396 + w1[u]] - 
   0.0361921 Log[0.95007 + w1[u]]) + 1.34588*10^-7 w1[u] == 8.74776*10^-6 
   - 1.34588*10^-7 u, w1[u]] *)

Thus, we must be satisfied with u as a function of w1.

su = u /. Simplify[Solve[%[[1]] /. w1[u] -> w1, u]][[1, 1]]
(* 64.9966 - 1. w1 - 1.04938 Log[-51.909 + w1] - 
   0.854463 Log[-30.0396 + w1] + 0.90384 Log[0.95007 + w1] *)
ParametricPlot[{su, w1}, {w1, 50, 72}, AxesLabel -> {u, w1}, AxesOrigin -> {0, 50}]

eqLRC can be solved in terms of ProductLog.