# How can I combine NDSolve[] with NSolve[]

Consider the two below ODEs with following dependent B.C. In the boundary conditions $k$ is an unknown.

$$y_1''(x)=x$$ $$y_2''(y)=0$$ $$y_1(1)=y_2(1)$$ $$y_1(-k)=y_2(-k)$$ $$\frac{\text{dy}_1(1)}{\text{dx}}=\frac{\text{dy}_2(1)}{\text{dx}}$$ $$\frac{\text{dy}_1(-k)}{\text{dx}}=\frac{\text{dy}_2(-k)}{\text{dx}}$$ $$\int_0^k {{y_2}\left( x \right)} dx = k$$ In this case, in order to find all unknowns, first, DSolve[] is used to solve the ODEs. Second, B.C. is applied to create 5 equations and unknowns. At the end, NSolve[] is implied to solve the system of equation with 5 unknowns. Following code is what I explained.

odey1 = DSolveValue[y1''[x] == x, y1[x], x];
odey2 = DSolveValue[y2''[y] == 0, y2[y], y] /. {C[1] -> C[3],
C[2] -> C[4]};
eqn1 = (odey1 /. x -> 1) == (odey2 /. y -> 1);
eqn2 = (odey1 /. x -> -k) == (odey2 /. y -> -k);
eqn3 = ((D[odey1, x]) /. x -> 1) == ((D[odey2, y]) /. y -> 1);
eqn4 = ((D[odey1, x]) /. x -> -k) == ((D[odey2, y]) /. y -> -k);
eqn5 = Integrate[odey2, {y, 0, k}] == k

coef = NSolve[
eqn1 && eqn2 && eqn3 && eqn4 && eqn5, {C[1], C[2], C[3], C[4], k}]
odey1 /. coef;
odey2 /. coef;


{{C[1] -> -7.95938*10^15 + 5.84722*10^15 I, C[2] -> -2.92835*10^15 + 1.97798*10^15 I, C[3] -> -7.95938*10^15 + 5.84722*10^15 I, C[4] -> -2.92835*10^15 + 1.97798*10^15 I, k -> -4.88246 + 1.65421 I}, {C[1] -> 8.65911*10^11 - 7.14483*10^11 I, C[2] -> 3.70956*10^11 - 2.62275*10^11 I, C[3] -> 8.65911*10^11 - 7.14483*10^11 I, C[4] -> 3.70956*10^11 - 2.62275*10^11 I, k -> -5.03678 + 0.40725 I}, {C[1] -> -2.25067*10^10 - 1.96804*10^10 I, C[2] -> 5.21094*10^9 - 1.92651*10^9 I, C[3] -> -2.25067*10^10 - 1.96804*10^10 I, C[4] -> 5.21094*10^9 - 1.92651*10^9 I, k -> 5.50713 + 9.1575 I}, {C[1] -> 4441.08 + 1806.41 I, C[2] -> 8880.03 + 3613.21 I, C[3] -> 4440.75 + 1806.41 I, C[4] -> 8880.53 + 3613.21 I, k -> -0.99988 - 2.69101*10^-6 I}}

My real case is much complicated and the ODEs doesn’t have close form. Therefore, I have to use NDSolve[] for solving ODEs. In this case I don’t know how can I combine NDSolve[] with NSolve[] to find the unknowns C[1], C[2], C[3], C[4] and $k$.

I will appreciate, if somebody can help me.

• You will probably have to use FindRoot. If only one solution is sought, it is probably the way to go. If you seek several solutions, then it might not; might depend on the actual problem. BTW, NDSolve[] implements a shooting method for BVPs -- would that work for your actual problem? – Michael E2 Jun 13 '16 at 10:49
• After looking more closely, I don't think the shooting method is set up to deal with such complicated BCs. -- You have some ells l in your code which I suspect are supposed to be ones 1. Is that right? – Michael E2 Jun 13 '16 at 11:51
• @MichaelE2 Thanks for your respond. I have modified the problem to have the solution. The problem that I presented before was a hypothetical example. Now, with current example do you have any Idea how can I change DSolve[] with NDSolve? Thanks. – Emad Jun 14 '16 at 20:19

## Preamble

After solving the differential equations it appears that there are five equations with five unknowns (four generated from the two differential equations and the parameter k.

However there are only four pieces of information in the five equations so there still is one arbitrary constant.

There are many routes that lead to this conclusion, perhaps the easiest to see is as follows. Since y''[x]=0 this implies that y'[x] is a constant. This means that the right hand side of eqns3 and eqns4 are identical.

As a result eqns3 and eqns4 collapse into one equation where the left hand sides are now equal.

## DSolve - boundary conditions

The equations will be re-written so that two boundary conditions (eqn1 and eqn3) will be incorporated into the equations. The motivation for this is to pave the way for a numerical solution.

The boundary conditions for the main function and its derivative are set at x = 1 and are equated to the arbitrary constants a and b.

odc1 = DSolveValue[{y1''[x] == x, y1[1] == a, y1'[1] == b}, y1[x], x]
(* 1/6 (2 + 6 a - 6 b - 3 x + 6 b x + x^3) *)

odc2 = DSolveValue[{y2''[x] == 0, y2[1] == a, y2'[1] == b}, y2[x], x]
(* a - b + b x *)


I find it useful to create functions for the solutions and (later) the derivatives.

yc1[x_, a_, b_] := (2 + 6 (a - b) + (6 b - 3 ) x + x^3)/6
yc2[x_, a_, b_] := a + b (x - 1)


Next take the derivative of yc1 and yc2.

D[yc1[x, a, b], x]
(* 1/6 (-3 + 6 b + 3 x^2) *)

D[yc2[x, a, b], x]
(* b *)


Observe that the derivative of yc2 is a constant. Again this is why two of the equations collapse.

Functions are defined for the derivatives

dyc1[x_, a_, b_] := (2 b - 1 + x^2)/2
dyc2[x_, a_, b_] := b


Next the other three equations are written. eqc1 is set to the condition where the two solutions are equated at x = -k.

eqc1 = yc1[-k, a, b] == yc2[-k, a, b]
(* 1/6 (2 + 6 (a - b) - (-3 + 6 b) k - k^3) == a + b (-1 - k) *)


eqc2 is set to the condition where the two derivates are equated at x = -k.

eqc2 = dyc1[-k, a, b] == dyc2[-k, a, b]
(* 1/2 (-1 + 2 b + k^2) == b *)


eqc3 is set to the condition where the integral is equated to k

Integrate[yc2[x, a, b], {x, 0, k}]
(* a k - b k + (b k^2)/2 *)

eqc3 = (a - b ) k + (b k^2)/2 == k;


Next solve these three equations allowing a to be an arbitrary constant

Solve[{eqc1, eqc2, eqc3}, {b, k}]
(* {{b -> 2/3 (-1 + a), k -> -1}} *)


The solution shows a relationship between a and b. k is found to be -1.

A second form for the functions yc1 and yc2 are created that honor the solution and only take the a parameter.

yc1[x_, a_] := Module[
{
b = (2 (a - 1))/3
},
(2 + 6 (a - b) + (6 b - 3 ) x + x^3)/6
]

yc2[x_, a_] := Module[
{
b = (2 (a - 1))/3
},
a + b (x - 1)
]


The function has the same name as the original but is distinct because of only two arguments.

Below is a plot for three values of a.

Plot[{yc1[x, -1], yc1[x, 0], yc1[x, 1]}, {x, -3, 3},
PlotStyle -> {Blue, Black, Red}]


Plot[{yc2[x, -1], yc2[x, 0], yc2[x, 1]}, {x, -3, 3},
PlotStyle -> {Blue, Black, Red}]


## Numerical Solution

We will attempt to do follow the same steps as above but use a numerical solver. In order to allow the differential equation to accept non-numerical parameters it is necessary to use ParametricNDSolve.

opc1 = ParametricNDSolveValue[{y1''[x] == x, y1[1] == a, y1'[1] == b},
y1, {x, -5, 5}, {a, b}]
opc2 = ParametricNDSolveValue[{y2''[x] == 0, y2[1] == a, y2'[1] == b},
y2, {x, -5, 5}, {a, b}]


The output is a ParametricFunction that we can use downstream. Functions are defined that require numerical inputs (again, needed later).

ynp1[x_?NumericQ, a_?NumericQ, b_?NumericQ] := opc1[a, b][x]
ynp2[x_?NumericQ, a_?NumericQ, b_?NumericQ] := opc2[a, b][x]


Here are the derivates

dynp1[x_?NumericQ, a_?NumericQ, b_?NumericQ] := opc1[ a, b]'[x]
ynp2[x_?NumericQ, a_?NumericQ, b_?NumericQ] := opc2[ a, b]'[x]


For the integral we will define a numerical form

int2[a_?NumericQ, b_?NumericQ, k_?NumericQ] :=
NIntegrate[ynp2[x, a, b], {x, 0, k}]


By trial and error I found that FindRoot wanted the a function as an input rather than lhs == rhs.

The functions are described by taking the difference between the left and right hand sides.

eqd1 = ynp1[-k, a, b] - ynp2[-k, a, b]
eqd2 = dynp1[-k, a, b] - dynp2[-k, a, b]
eqd3 = int2[a, b, k] - k


FindMinimum worked by asking it to minimize the sum of the squares of the above expressions but the parameters a was set to an arbitrary constant (unity for this example). Surprisingly FindMinimum would not accept eqd1 through eqd3, one had to manually insert the corresponding expression.

sol = FindMinimum[
(ynp1[-k, 1, b] - ynp2[-k, 1, b])^2 +
(dynp1[-k, 1, b] - dynp2[-k, 1, b])^2 +
(int2[1, b, k] - k)^2,
{{b, -1}, {k, -1}}
]

(* {1.97215*10^-31, {b -> 9.73186*10^-16, k -> -1.}} *)


FindRoot was a bit fussy in that it complained if the number of equations and unknowns were not identical. eqd1 was replaced with the statement substracting a from the desired value. Interesting that FindRoot happily accepted eqd2 and eqd3.

FindRoot[
{
a - 1,
eqd2,
eqd3
},
{{a, 1}, {b, 0}, {k, -1}}
]
(* {a -> 1., b -> 0., k -> -1.} *)


Both of these results match the DSolve.

A function to plot ynp1 with the solution constants can be defined as:

plotYp1[astart_?NumericQ, color_] := Module[
{
sol = FindMinimum[
(ynp1[-k, astart, b] - ynp2[-k, astart, b])^2 +
(dynp1[-k, astart, b] - dynp2[-k, astart, b])^2 +
(int2[astart, b, k] - k)^2,
{{b, -1}, {k, -1}}
]
},

Plot[(ynp1[x, astart, b] /. sol[[2]]), {x, -3, 3},
PlotRange -> {{-3.1, 3.1}, {-2, 4.65}},
PlotStyle -> color]
]


And a duplicate of the DSolve plot is made.

Show[
plotYp1[-1, Blue],
plotYp1[0, Black],
plotYp1[1, Red]
]


Similar results apply for the second function:

plotYp2[astart_?NumericQ, color_] := Module[
{
sol = FindMinimum[
(ynp1[-k, astart, b] - ynp2[-k, astart, b])^2 +
(dynp1[-k, astart, b] - dynp2[-k, astart, b])^2 +
(int2[astart, b, k] - k)^2,
{{b, -1}, {k, -1}}
]
},

Plot[(ynp2[x, astart, b] /. sol[[2]]), {x, -3, 3},
PlotRange -> {{-3.1, 3.1}, {-2, 4.65}},
PlotStyle -> color]
]

Show[
plotYp2[-1, Blue],
plotYp2[0, Black],
plotYp2[1, Red]
]


Hope this helps you with your actual problem. Good luck!

• I really appreciate you. Your answer is so professional. I am trying to digest what you did and implement your procedure to my actual problem. Thanks a lot. – Emad Jun 17 '16 at 3:15