# Solve BVP involving user defined function [closed]

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?

• – Michael E2 Mar 20 '18 at 18:24
• 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$. – Michael E2 Mar 20 '18 at 18:27
• Did you try using NDSolve, and if so, what did you do and what went wrong? – Carl Woll Mar 20 '18 at 19:09
• 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 – Artes Mar 21 '18 at 10:02
• Moreover take a look at this answer. – Artes Mar 21 '18 at 10:03

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 == c, y[L] == 0}, y, {x, 0, L}]

Plot[Evaluate[y[x] /. ndsol], {x, 0, 10}] Show the error of the solution

Plot[Evaluate[y''[x] + f[y[x]] /. ndsol], {x, 0, 10},
PlotRange -> 10^-4] 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)]


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], {K, 1, K}] + C],
{K, 1, y[x]}]^2 == (C + 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], {K, 1, K}] + C],
{K, 1, y[x]}] == C + x /. {C -> c1, C -> c2, K -> k1,
K -> 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 -> 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].