# NDSolving differential equations with complicated initial conditions

I want to solve numerically the following system of differential equations involving three functions $f[x],g[x],h[x]$ where $x\in \mathbb{R}$

$$f''[x] -4 g'[x]\,f'[x]=0\\ g''[x]-\lambda\, h'[x]^2=0\\ g'[x]^2 f[x] -\frac{1}{4}\, f'[x]\, g'[x]-1-\lambda\, h'[x]^2=0$$

where $\lambda \ll 1$ and the variable $x\in [0,1]$. I want to impose the following initial conditions

$$f[0]=1\,,\quad f[1]=0\,,\quad g[0]=0\,,\quad g[1]=-\ln\left(\frac{4e^{-1}}{|\,f'[1]\,|}\right)\,,\quad h[0]=1.$$

This system can be solved analytically for $\lambda=0$ and the solution is determined up to $h[x]$

$$g[x] = x\,,\qquad f[x] = 1-e^{4(x-1)}\,.$$

I want to generate the solutions for a finite set of values of $\lambda \ll 1$. However, using NDSolve, I don't know how to implement the boundary condition on $g[1]$ which depends on the derivative of $f[x]$ computed in $x=1$. My attempt for $\lambda=0.1$ (which obviously fails) is the following

\[Lambda] = 0.1;
NDSolve[{f''[x] - 4 g'[x] f'[x] == 0, g''[x] - \[Lambda] h'[x]^2 == 0, g'[x]^2 f[x] - 1/4 f'[x] g'[x] - 1 - \[Lambda] h'[x]^2 == 0, f[0] == 1, f[1] == 0, g[0] == 0, h[0] == 1, g[1] == -Log[4 E^-1/Abs[f'[1]]]}, {f[x], g[x], h[x]}, {x, 0, 1}]


I get the following error

I understand the error is probably related to the fact that the boundary condition on $g[1]$ is written in a wrong way. Indeed, $f'[1]$ is not even define before to do the numerical computation.

Questions:

• How can I solve the system for $\lambda\ll 1$ but still $\lambda \neq 0$?
• How can I extend the code to compute the solution for $\lambda=0$? Indeed, in this case the system is "apparently" overdetermined because we have three equations and only two variables ($f[x]$ and $g[x]$). However, the first equation is redundant because it can be derived by differentiation of the third equation. If I set $\lambda=0$ in my code, I get the error
• Try simplifying the equations by eliminating h and lamda between the second and third equations. Now you have two equations in two unknowns, which should be easier to solve. In addition, invert your fourth boundary condition, so that it expresses f' in terms of g. By the way, your analytical solution for lamda == 0 does not satisfy f[0] == 1. – bbgodfrey Sep 11 '18 at 4:55

This ODE system can be solved as outlined in my earlier comment above. For convenience, define

eq1 = f''[x] - 4 g'[x] f'[x];
eq2 = g''[x] - λ h'[x]^2;
eq3 = g'[x]^2 f[x] - 1/4 f'[x] g'[x] - 1 - λ h'[x]^2;


First, eliminate λ h'[x]^2 between eq2 and eq3

eq4 = Simplify[eq2 - eq3]
(* 1 + 1/4 f'[x] g'[x] - f[x] g'[x]^2 + g''[x] *)


Next, solve the fourth boundary condition for f'[1]. The result, of course, depends on the sign of f'[1]. If positive,

nbc0 = Solve[g[1] == -Log[4 E^-1/f'[1]], f'[1], Reals][[1, 1]] /. Rule -> Equal
(* f'[1] == 4 E^(-1 + g[1]) *)


and if negative,

nbc2 = Solve[g[1] == -Log[-4 E^-1/f'[1]], f'[1], Reals][[1, 1]] /. Rule -> Equal
(* f'[1] == -4 E^(-1 + g[1]) *)


NDSolve then solves the ODEs without difficulty. For positive f'[1],

s0 = NDSolve[{eq1 == 0, eq4 == 0, f[0] == 1, f[1] == 0, g[0] == 0, nbc0},
{f[x], g[x]}, {x, 0, 1}] // Values // Flatten;

Plot[First[s0], {x, 0, 1}, PlotRange -> All, AxesLabel -> {x, f},
LabelStyle -> {Bold, Black, Medium}, ImageSize -> Large]
Plot[Last[s0], {x, 0, 1}, PlotRange -> All, AxesLabel -> {x, g},
LabelStyle -> {Bold, Black, Medium}, ImageSize -> Large]


The parallel computation for negative f'[1] (i.e., with nbc2) yields

Although h[x] now can be determined for specified values of λ in a straightforward manner by applying NDSolve to eq2, it is more instructive to write the solution of eq2 as

s1 = 1 + λ^-(1/2) NDSolveValue[{h'[x] == Evaluate@Sqrt[D[Last[s0], {x, 2}]],
h[0] == 0}, h[x], {x, 0, 1}]


(A second solution to eq2 is 1 - λ^-(1/2) ….) Visibly, h[x] becomes infinite as λ approaches 0. Moreover, numerical examination of g''[x]] shows that it negative in all cases, causing h[x] to be purely imaginary. A plot of the integral in s1 for positive f'[1] is

hint = s1[[2, 2]]
Plot[Im@hint, {x, 0, 1}, PlotRange -> All, AxesLabel -> {x, "hint"},
LabelStyle -> {Bold, Black, Medium}, ImageSize -> Large]


and similarly for negative f'[1].

Addendum - Symbolic Solution in Question

Also as noted in my earlier comment, the symbolic solution for λ == 0 given in the question is incorrect, because it gives f[0] == 1 - 1/E^4, not 1 as required. Interestingly, if f[0] == 1 - 1/E^4 were the desired boundary condition, then the entire problem could be solved symbolically for negative f'[1] with f[x] and g[x] as in the question and h[x] == 1 for all values of λ.

This method is a little more involved, but it does not require eliminating one of the equations, which may not always be an option with a slightly different problem.

Clear["Global*"]


Solve the first equation to get f in terms of g.

ode1 = f''[x] - 4 g'[x] f'[x] == 0

DSolve[ode1, f[x], x] // Flatten

f[x_] = f[x] /. % /. {C[1] -> c1, C[2] -> c2}

(*Integrate[c1*E^(4*g[K[1]]), {K[1], 1, x}] + c2*)


Apply a bc

f[1] == 0
(*c2 == 0*)

c2 = 0

f[x]
(*Integrate[c1*E^(4*g[K[1]]), {K[1], 1, x}]*)


I choose a more readable form and pull c1 out of the integral.

f[x_] = c1*Integrate[E^(4*g[z]), {z, 1, x}];


Apply the other bc on f

f[0] == 1;


and solve for c1.

Solve[%, c1] // Flatten
(*{c1 -> 1/Integrate[E^(4*g[z]), {z, 1, 0}]}*)

Clear[g]


The condition for g[1]

eq = g[1] == -Log[4/(E*Abs[Derivative[1][f][1]])]
(*g[1] == -Log[(4*E^(-(4*Re[g[1]]) - 1))/Abs[c1]]*)

\$Assumptions = c1 \[Element] Reals && g[1] \[Element] Reals

Solve[eq, g[1]] // Flatten // Simplify
(*{g[1] -> ConditionalExpression[1/3 (Log[4/Abs[c1]] - 1),
C[1] \[Element] Integers && C[1] == 0]}*)

(g[1] /. %[[1]])[[1]]
(*1/3 (Log[4/Abs[c1]] - 1)*)


Use the above for the condition for g[1] and iterate for c1.

g1 := 1/3 (Log[4/Abs[c1]] - 1)

\[Lambda] = .1;


Starting value for c1

c1 = 1;

Do[
Clear[f, g, h];
c2 = c1;
sol = NDSolve[{Derivative[2][f][x] - 4*Derivative[1][g][x]*Derivative[1][f][x] == 0,
Derivative[2][g][x] - \[Lambda]*Derivative[1][h][x]^2 == 0,
Derivative[1][g][x]^2*f[x] - (1/4)*Derivative[1][f][x]*Derivative[1][g][x] - 1 -
\[Lambda]*Derivative[1][h][x]^2 == 0, f[0] == 1, f[1] == 0, g[0] == 0, h[0] == 1,
g[1] == (1/3)*(Log[4/Abs[c1]] - 1)}, {f[x], g[x], h[x]}, {x, 0, 1}];

f[x_] = f[x] /. sol[[1]];
g[x_] = g[x] /. sol[[1]];
h[x_] = h[x] /. sol[[1]];
int = N[Integrate[E^(4*g[x]), {x, 1, 0}]];
c1 = 1/int;
If[Abs[c2 - c1] < 10^(-5), Break[]];
, 500]


The above iterates until the change in c1 is < 10^(-5)

Plots

Plot[f[x], {x, 0, 1}]


Plot[g[x], {x, 0, 1}]


h is complex everywhere except at x = 0. h'[x]^2 is real and probably could be replaced in the equations by a different function of x. You would need to know the bc for that function, however.

Check the bc's

f[0] // Chop
(*1.*)

f[1] // Chop
(*-4.80252*10^-7*)

g[0] // Chop
(*0*)

g[1] + Log[4/(E Abs[f'[1]])] // Chop
(*0.0000986539*)

h[0] // Chop
(*1.*)


which is the only x value for which h is real.

Check the differential equations

Plot[f''[x] - 4 g'[x] f'[x], {x, 0, 1}, PlotRange -> {-.1, .1}]


Plot[g''[x] - \[Lambda] h'[x]^2, {x, 0, 1}, PlotRange -> {-.1, .1}]


Plot[g'[x]^2 f[x] - 1/4 f'[x] g'[x] - 1 - \[Lambda] h'[x]^2, {x, 0,
1}, PlotRange -> {-.1, .1}]


It is interesting in this case that f and g do not depend on lambda at all. The case lambda = 0` cannot be solved in this case because the system is over determined. You can use any two equations to solve for f and g, but they will not satisfy the unused third equation.

This particular iterative solution goes pretty fast, but for more complicated cases, it could be very slow because the entire system needs to be solved in each iteration.