4
$\begingroup$

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!

Improve this question
$\endgroup$
1
  • $\begingroup$ Thank you @Xavier! I updated the typos in my post. $\endgroup$
    – user38307
    Commented Sep 11, 2016 at 1:56

1 Answer 1

Reset to default
3
$\begingroup$

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

enter image description here

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]

enter image description here

The overlap is smooth, as desired.

Addendum

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.

Improve this answer
$\endgroup$
8
  • $\begingroup$ 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? $\endgroup$
    – user38307
    Commented Sep 11, 2016 at 20:32
  • $\begingroup$ @user38307 I copied the code from my answer into a new notebook and ran it with Mathematica 10.4.1 and with 11.0, obtaining the same results as before. What version of Mathematica are you using? $\endgroup$
    – bbgodfrey
    Commented Sep 11, 2016 at 23:27
  • $\begingroup$ @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. $\endgroup$
    – bbgodfrey
    Commented Sep 12, 2016 at 15:59
  • $\begingroup$ Upgrading to version 11 reproduces your solution just great! I'm trying to plot the "full" solution (i.e. from x1 to x2) and I'm not having much luck. In trying to understand your approach, I think s1 solves the 2 Diff eqs separately using the "known" boundary conditions leaving 2 unresolved constants for each sol. s2 solves the matching conditions w'=w2' at xg. s3 appears to evaluate xg s.t. s3==0 gives xg. Do you not need the higher order matching condition (w''')? It would be great if you could include the way to get the full solution to w, w2 in your excellent answer. $\endgroup$
    – user38307
    Commented Sep 16, 2016 at 16:49
  • $\begingroup$ @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. $\endgroup$
    – bbgodfrey
    Commented Sep 16, 2016 at 16:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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