3
$\begingroup$

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.

|improve this question
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Alexei Boulbitch Nov 4 '15 at 13:19
  • $\begingroup$ @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." $\endgroup$ – jak Nov 4 '15 at 13:24
  • $\begingroup$ @AlexeiBoulbitch I updated the question with my complete problem, as the previously stated part of the problem was unsolvable $\endgroup$ – jak Nov 4 '15 at 13:38
  • $\begingroup$ Now it seems to be one extra boundary condition. Indeed, 6 equations and 7 boundary conditions. $\endgroup$ – Alexei Boulbitch Nov 4 '15 at 13:52
  • $\begingroup$ @AlexeiBoulbitch Again thanks. I removed one of them. $\endgroup$ – jak Nov 4 '15 at 13:56
2
$\begingroup$

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}]

enter image description here

eqLRC can be solved in terms of ProductLog.

|improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.