0
$\begingroup$

I have a second order ODE which I would like to solve numerically:

$\frac{d^2 y}{dx^2}+f(y)=0$

BCs: $ y(0) = c$, $y(L)=0$

Where $c$ is a known constant and $f(y)$ is a complicated user defined function, which includes the use of FindRoot.

Naturally, this can be written as a first order system:

$\frac{dy}{dx}=\rho$

$\rho\frac{d\rho}{dy} + f(y) = 0$

If $f(y)$ was a nice simple function this could be solved fairly easily since it is separable. However, I'm not sure how to proceed in the case where $f(y)$ is a user defined function which is not known in closed form.

Any suggestions of how to proceed?

Improve this question
$\endgroup$
5
  • 2
    $\begingroup$ See mathematica.stackexchange.com/questions/163664/… $\endgroup$
    – Michael E2
    Commented Mar 20, 2018 at 18:24
  • 1
    $\begingroup$ If the answers there do not work in you particular case, a working example for f[y] would make it easier for people to try out their ideas. Also, since it is separable, one can use NDSolve to integrate $\int f(y)\; dy$. $\endgroup$
    – Michael E2
    Commented Mar 20, 2018 at 18:27
  • 2
    $\begingroup$ Did you try using NDSolve, and if so, what did you do and what went wrong? $\endgroup$
    – Carl Woll
    Commented Mar 20, 2018 at 19:09
  • $\begingroup$ In case of $f(y)=c_1 y$ you have a simple solution in terms of trygonometric functions, for $f(y)= a_3 y^3+a_2 y^2+a_1 y+ a_0$ this can be solved in terms of elliptic functions, for example see this or this answer $\endgroup$
    – Artes
    Commented Mar 21, 2018 at 10:02
  • $\begingroup$ Moreover take a look at this answer. $\endgroup$
    – Artes
    Commented Mar 21, 2018 at 10:03

3 Answers 3

Reset to default
2
$\begingroup$

As an example suppose f[y] is the root of z^3 + 20 Cos[y] - 1/2 

f[y_?NumericQ] := 
        z /. First@FindRoot[z^3 + 20 Cos[y] - 1/2 == 0, {z, 1}]

L = 10; ndsol[c_] := 
          NDSolve[{y''[x] + f[y[x]] == 0, y[0] == c, y[L] == 0}, y, {x, 0, L}]

Plot[Evaluate[y[x] /. ndsol[1]], {x, 0, 10}]

enter image description here

Show the error of the solution

Plot[Evaluate[y''[x] + f[y[x]] /. ndsol[1]], {x, 0, 10}, 
    PlotRange -> 10^-4]

enter image description here

Improve this answer
$\endgroup$
0
$\begingroup$

Since you have already defined the function as f[y], this implies that it can be expressed in terms of y[x]. So, here is what I suggest:

  1. After defining your reference coordinate system, and your domain x-y, tabulate (and plot) f[y] as a function of y[x].

  2. Assuming the function is single-valued, do a curve fit of f[y] versus y[x] to obtain the closed form expression for f[y]. If f[y] is multivalued, you may need to parametrize f[y] in terms of a new parameter s[x] or t[y], before you proceed. If you are unable to parametrize the multivalued function f[y], you can divide the domain into sub-domains where the function is single-valued before you proceed to solve the ODE separately for each sub-domain.

  3. Assuming that f[y] is single-valued along the path y[x] in the domain x-y, you can obtain the best fit expression for example, f[y] = a1 y[x] + a2 y[x]^2 + a3 y[x]^3 + ... Please note that a polynomial equation may not be the optimal expression for your fit. You need to use some judgement here, using the plot as a guide.

  4. Now you are ready to solve the original second order ODE directly, using DSolve, DSolveValue or NDSolve, and the two given boundary conditions.

Example: For the simple case where f[y] is a linear function of y[x], the best curve fit yields f[y]=a1*y[x], and given the two boundary conditions, the solution becomes:

y[x] = c*Csc[L]*Sin[(L-x)]
Improve this answer
$\endgroup$
0
$\begingroup$

We can at least partially solve this if f[y] is integrable.

ode = y''[x] + f[y[x]] ==0

DSolve[ode, y[x], x]

(* Solve[Integrate[1/Sqrt[2*Integrate[-f[K[1]], {K[1], 1, K[2]}] + C[1]], 
     {K[2], 1, y[x]}]^2 == (C[2] + x)^2, y[x]] *)

Look at the equation to solve and simplify a few things. Take the positive Sqrt of each side as a starting point.

eq = Integrate[1/Sqrt[2*Integrate[-f[K[1]], {K[1], 1, K[2]}] + C[1]], 
     {K[2], 1, y[x]}] == C[2] + x /. {C[1] -> c1, C[2] -> c2, K[1] -> k1, 
    K[2] -> k2}
(* Integrate[1/Sqrt[c1 + 2*Integrate[-f[k1], {k1, 1, k2}]], 
   {k2, 1, y[x]}] == c2 + x *)

Apply the bc's

eq /. x -> 0 /. y[0] -> c
(* Integrate[1/Sqrt[c1 + 2*Integrate[-f[k1], {k1, 1, k2}]], 
   {k2, 1, c}] == c2 *)

eq /. x -> L /. y[L] -> 0
(* Integrate[1/Sqrt[c1 + 2*Integrate[-f[k1], {k1, 1, k2}]], 
   {k2, 1, 0}] == c2 + L *)

We could possibly solve these equations for c1 and c2 if we can integrate f[y] either algebraically or numerically, but it is hard to go any further without knowing f[y].

Improve this answer
$\endgroup$

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