Solve Differential Equations with Boundary Condition at Variable Point

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.

• You seem to need one more boundary condition, since w2R[M]==w2L[M]==w4C[M]  is equivalent to w2R[M]==w2L[M]  and w2L[M]==w4C[M] , while you have three dependent variables. Commented Nov 4, 2015 at 13:19
• @AlexeiBoulbitch Of course, thanks! I have another function $w1$ (from a solution of NDsolve] and the other boundary condition is that at some other variable point $M2$ I should have $w1[M2]=w2R[M2]+w4C[M2]$. Thus the solution should depend on the two parameters $M2$ and $M$. Unfortunately, Mathematica then tells me: ""There are fewer dependent variables, {w2L[u],w2R[u],w4C[u]}, than \ equations, so the system is overdetermined."
– jak
Commented Nov 4, 2015 at 13:24
• @AlexeiBoulbitch I updated the question with my complete problem, as the previously stated part of the problem was unsolvable
– jak
Commented Nov 4, 2015 at 13:38
• Now it seems to be one extra boundary condition. Indeed, 6 equations and 7 boundary conditions. Commented Nov 4, 2015 at 13:52
• @AlexeiBoulbitch Again thanks. I removed one of them.
– jak
Commented Nov 4, 2015 at 13:56

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.

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.