# Solving coupled Differential equations with matching condition

I am somewhat new to using Mathematica and was wondering if it can solve the following: $(a \rightarrow h)$ are constants

$a w^4[x] =-b + c(d-w[x]-e x),\hspace{2em} x_1 \leq x \leq x_g \\ a w^4[x] =-b + f(d-w[x]),\hspace{2em} x_g < x \leq x_2$

With the boundary conditions:

$w[x_1]=d - g + e x_1\\ w'[x_1]=e \\ [w]=[w']=[w'']=[w''']=0, \; x=x_g \; \textrm{ (jump conditions)} \\ w[x_g]=d+e x_g \\ w''[x_2]=w'''[x_2]=0$

There is also a constraint: $\int_{x1}^{x2} \sqrt{1+w'^2} dx=h$

I believe this model is analogous to solving beam equations with matched conditions except it is nonlinear: basically a cantilevered beam partially resting on a foundation with a pinning point at the end of the foundation that is allowed to move.

I tried searching the site for using jump conditions, but I didn't have much luck with the answers from the top hits. I was thinking that I could decouple the equations at $x_g$ and include matching conditions based on 2 differential equations: $w_1[x], w_2[x]$. However, the length constraint doesn't appear to be applicable.

So, I tried the following:

a = 8.06*10^16;
b = 8.99*10^6;
c = 8.99*10^11;
d = 500.;
e = -.268;
f = 10074.4;
g = 10.^(-5);
x1 = 0.;
x2 = 10000.;

NDSolve[{
a*w''''[x] == -b + c (d - w[x] + e x),
a*w2''''[x] == -b + f (d - w2[x]),
w[x1] == d - g + e x1,
w'[x1] == e,
w[xg] == w2[xg],
w'[xg] == w2'[xg],
w''[xg] == w2''[xg],
w'''[xg] == w2'''[xg],
w2[xg] == d + e xg,
w2''[x2] == 0,
w2'''[x2] == 0
},
{w, w2, xg}, {x, x1, x2}
]


I get the error:

NDSolve::ndsv: Cannot find starting value for the variable w''


I think this error is pointing to a missing boundary condition. If I specify xg, I get another error, which I'm not sure how to interpret. I would greatly appreciate any comments on solving these equations. Thanks!

• Thank you @Xavier! I updated the typos in my post. – user38307 Sep 11 '16 at 1:56

This is an eigenvalue problem of sorts, with xg the quantity to be determined. The strategy is to determine w and w2 using DSolve and all the boundary conditions except w2[xg] == d + e xg, which is employed only at the end of the computation to determine xg. To avoid loss-of-precision issues, all constants must be rationalized,

rules = Rationalize[{a -> 8.06*10^16, b -> 8.99*10^6, c -> 8.99*10^11,
d -> 500, e -> -.268, f -> 10074.4, g -> 10^-14}, 10^-5];


and even then none of the constants except x1 and x2 given their numerical values at first. Unfortunately, with all boundary conditions (except w2[xg] == d + e xg) included, DSolve crashes the kernel after running a long time. So, begin by solving

s1 = Collect[DSolve[{a*w''''[x] == -b + c (d - w[x] + e x),
a*w2''''[x] == -b + f (d - w2[x]), w[x1] == d - g + e x1,
w'[x1] == e, w2''[x2] == 0, w2'''[x2] == 0} /. {x1 -> 0, x2 -> 10000},
{w[x], w2[x]}, x] // Flatten, {C[1], C[2], C[5], C[6]}, Simplify];


which produces an expression a bit long to be reproduced here. Next, use the four remaining boundary conditions to determine the four constants of integration, {C[1], C[2], C[5], C[6]}

s2 = Solve[Collect[D[(w[x] == w2[x]) /. s1, {x, #}] /. x -> xg, {C[1], C[2], C[5], C[6]},
Simplify] & /@ Range[0, 3], {C[1], C[2], C[5], C[6]}] // Flatten;


which is enormous. Finally, the expression determining xg is

s3 = (w2[x] - (d + e xg)) /. s1 /. x -> xg /. s2 /. rules;


A plot of s3 give the approximate value of xg.

ListLinePlot[Table[Chop[N[s3, 30]], {xg, 950, 1010}], DataRange -> {950, 1010}]


Without a precision of 30, N[s3] is highly inaccurate near xg == 980. Finally, xg can be computed by

FindRoot[Chop[N[s3, 30]], {xg, 976, 980}, WorkingPrecision -> 30, Evaluated -> False]
(* {xg -> 977.546450674047397721034539858} *)


Alternative (Better) Approach

In the computation above, DSolve gives an expression containing complex numbers. Although the imaginary parts eventually cancel, as they should, higher precision along with Chop is needed to achieve the cancellation. This alternative approach transforms the DSolve expression to eliminate the complex numbers at the outset. (It does not seem possible, however, to instruct DSolve itself not to include explicitly complex numbers in its output. See question 126132 and associated comments.)

Collect[DSolve[{a*w''''[x] == -b + c (d - w[x] + e x),
a*w2''''[x] == -b + f (d - w2[x]), w[x1] == d - g + e x1,
w'[x1] == e, w2''[x2] == 0, w2'''[x2] == 0}, {w[x], w2[x]}, x,
Assumptions -> a > 0 && c > 0 && f > 0] // Flatten, { C[1], C[2],
C[5], C[6]}, Simplify];
Solve[{c1 == C[1] + C[2], c2 == I (C[1] - C[2]), c5 == C[5] + C[6],
c6 == I (C[5] - C[6])}, {C[1], C[2], C[5], C[6]}] // Flatten;
s4 = Collect[ComplexExpand[%% /. %, TargetFunctions -> {Re, Im}] /.
x z_ :> x FullSimplify[z, a > 0 && c > 0 && f > 0] /.
x1 z_ :> x1 FullSimplify[z, a > 0 && c > 0] /.
x2 z_ :> x2 FullSimplify[z, a > 0 && f > 0], {c1, c2, c5, c6},
Simplify] /. {x1 -> 0, x2 -> 10000}

(* {w[x] -> -c2 E^(-(((c/a)^(1/4) x)/Sqrt[2])) (-1 + E^(Sqrt[2] (c/a)^(1/4) x))
Sin[((c/a)^(1/4) x)/Sqrt[2]] + c1 E^(-(((c/a)^(1/4) x)/Sqrt[2]))
((1 - E^(Sqrt[2] (c/a)^(1/4) x)) Cos[((c/a)^(1/4) x)/Sqrt[2]] + 2 E^(Sqrt[2]
(c/a)^(1/4) x) Sin[((c/a)^(1/4) x)/Sqrt[2]]) + (-b + c d + c e x +
E^(((c/a)^(1/4) x)/Sqrt[2]) (b - c g) Cos[((c/a)^(1/4) x)/Sqrt[2]] -
E^(((c/a)^(1/4) x)/Sqrt[2]) (b - c g) Sin[((c/a)^(1/4) x)/Sqrt[2]])/c,
w2[x] -> d - b/f + c5 E^(-(((f/a)^(1/4) (20000 + x))/Sqrt[2]))
((E^(10000 Sqrt[2] (f/a)^(1/4)) + E^(Sqrt[2] (f/a)^(1/4) x))
Cos[((f/a)^(1/4) x)/Sqrt[2]] + E^(Sqrt[2] (f/a)^(1/4) x) (Sin[((f/a)^(1/4)
(-20000 + x))/Sqrt[2]] - Sin[((f/a)^(1/4) x)/Sqrt[2]])) +
c6 E^(-(((f/a)^(1/4) (20000 + x))/Sqrt[2])) (E^(Sqrt[2] (f/a)^(1/4) x)
Cos[((f/a)^(1/4) (-20000 + x))/Sqrt[2]] + E^(Sqrt[2] (f/a)^(1/4) x)
Cos[((f/a)^(1/4) x)/Sqrt[2]] + (E^(10000 Sqrt[2] (f/a)^(1/4)) + E^(
Sqrt[2] (f/a)^(1/4) x)) Sin[((f/a)^(1/4) x)/Sqrt[2]])} *)


Now, proceeding as above,

s5 = Solve[Collect[D[(w[x] == w2[x]) /. s4, {x, #}] /. x -> xg, {c1, c2, c5, c6},
Simplify] & /@ Range[0, 3], { c1, c2, c5, c6}] // Flatten;
s6 = (w2[x] - (d + e xg)) /. s4 /. x -> xg /. s5 /. rules;
ListLinePlot[Table[N[s6, 15], {xg, 950, 1010}], DataRange -> {950, 1010}]


yields the same curve as above but, as promised with less precision required and Chop not employed. Nearly the same result is obtained for xg.

(* {xg -> 977.546444338072792851682553074} *)


The total solution is easily plotted.

p = ListLinePlot[Table[N[w7], {x, 0, 1200, 10}], DataRange -> {0, 1200}];
p2 = ListLinePlot[Table[N[w27], {x, 500, 10000, 100}], DataRange -> {500, 10000}];
Show[p, p2, PlotRange -> All]


The overlap is smooth, as desired.

The parameter g in the question was changed by the OP from 10^-14 to 10^-5. No significant change in the results occurs. In fact, no significant change in the results occurs for g as large as 10.
• Thank you for this analysis @bbgodfrey. I copied the posted code, but I'm not getting a solution for s2... just {} . Do I need to set something in the preferences first? – user38307 Sep 11 '16 at 20:32
• @user38307 Although I no longer have a copy of 10.2, I am surprised that it would not compute s2. However, it is quite possible that it miscalculated s1, and the s2 calculation then would have failed too. In any case it makes sense to upgrade. There have been many improvements to DSolve (and other functions) recently. – bbgodfrey Sep 12 '16 at 15:59
• @user38307 As x moves toward x2, ever higher WorkingPrecision is needed. The s2 calculation uses all four matching conditions to determine the remaining four constants of integration. I do not understand your final sentence. By the way, I think that I have a slightly different approach to solving the equations that does not require such high WorkingPrecision. – bbgodfrey Sep 16 '16 at 16:55