I have a nonlinear boundary ODE,

$$y'' + 3 y y' = 0, \qquad y(0) = 0, y(2)=1 $$

I want to solve this using the finite-difference method. I obtained the result for the data set of $y(x)$ (including the boundary points),

$$(y_0,y_1,y_2,y_3,y_4,y_5,y_6,y_7,y_8,y_9,y_{10}) = (0, 0.318993138, 0.581636081, 0.7641886486, 0.875871387, 0.938497028, 0.9714644087, 0.9879847664, 0.9960488292, 1)$$

I know I could obtain the approximate function for $y(x)$ through interpolation. I'm interested in obtaining the derivative function $y'(x)$, I know I could just take the derivative of the interpolation function of $y(x)$, but the error is not small. So I'm thinking maybe there is a way to get the data points of $y'(x)$ from $y(x)$ because that might be more precise, but I'm not sure how. In textbooks, there are no formulas or instructions on how to obtain the data points of the derivative $y'$, i.e. $(y'_0,y'_1,y'_2,y'_3,y'_4,y'_5,y'_6,y'_7,y'_8,y'_9,y'_{10})$

The code is written in Mathematica.

(*Finite-difference method*)
(*a,b are endpoints, n is the grid points number, h is the step size*)
a = 0;
b = 2;
n = 9;
h = (b - a)/n;
alpha = 0;
beta = 1;
y[0] = alpha;
y[n] = beta;
eq = y''[x] + 3 y[x] y'[x];
For[i = 1, i <= n - 1, i++, eqn[i] = Simplify[Collect[eq /. {y''[x] -> ((y[i + 1] - 2 y[i] + y[i - 1])/h^2), y'[x] -> ((y[i + 1] - y[i - 1])/(2 h)), y[x] -> (y[i + 1] + y[i - 1])/2}, y[i]]]; Print["eqn[", i, "] = ", eqn[I]]];

(*Solving the nonlinear system of equations through Newton's method*)
(*Fx is the system of equations vector evaluated at xr[j], DFx is the Jacobian of the system of equations evaluated at xr[j]*)
m = 3;
For[j = 0, j <= m, j++, x[0] = {0, 0, 0, 0, 0, 0, 0, 0}; xr[j] = MapThread[#1 -> #2 &, {Array[y, Length[Table[i, {i, 1, n - 1}]]], x[j]}]; DFx = ResourceFunction["JacobianMatrix"][Table[eqn[i], {i, 1, n - 1}], Table[y[i], {i, 1, n - 1}]] /. xr[j]; Fx = Table[{eqn[i]}, {i, 1, n - 1}] /. xr[j]; x[j + 1] = x[j] + LinearSolve[DFx, -Fx] // Flatten]
datapoints = x[m] // N;
data = Join[{alpha}, datapoints, {beta}] // N

(*Compare Plots*)

(*Approximate solution curve y[x] through interpolation*)
a = 0;
b = 2;
n = 9;
h = (b - a)/n;
data = {0., 0.318993138, 0.581636081, 0.7641886486, 0.875871387, 0.938497028, 0.9714644087, 0.9879847664, 0.9960488292, 1.};
yfunc = Interpolation[MapThread[{#1, #2} &, {Table[i, {i, 0, 2, h}], data}], InterpolationOrder -> 3, Method -> "Spline"]

Plot[{yfunc[x], yfunc'[x]}, {x, 0, 2}, Frame -> True]

(*Solution curve y[x] computed internally, more precise?*)
ys[c_?NumericQ] := 
 Module[{e = c, s1}, 
  s1 = NDSolveValue[{v'[x] == (-3 y[x] y'[x]) /. y'[x] -> v[x], 
     y'[x] == v[x], y[0] == 0, y'[0] == e}, y[2], {x, 0, 2}]; s1]
sol = FindRoot[ys[c] == 1, {c, 2}, WorkingPrecision -> 10];
sols = NDSolveValue[{v'[x] == (-3 y[x] y'[x]) /. y'[x] -> v[x], 
    y'[x] == v[x], y[0] == 0, y'[0] == c /. sol}, y, {x, 0, 2}];

Plot[{Evaluate[y[x] /. y -> sols], Evaluate[y'[x] /. y -> sols]}, {x, 
  0, 2}, Frame -> True, PlotRange -> All]

(*Error Visualization*)
Plot[{Evaluate[y[x] /. y -> sols] - yfunc[x]}, {x, 0, 2}, 
 Frame -> True, PlotRange -> All]
Plot[{Evaluate[y'[x] /. y -> sols] - yfunc'[x]}, {x, 0, 2}, 
 Frame -> True, PlotRange -> All]

Plot of the solution curve from finite-difference through interpolation Image1

Plot of the solution curve from internal calculation of Mathematica using NDSolve Image2

$y(x)$ difference of the two approach


$y'(x)$ difference of the two approach


As you see the error of $y'(x)$ is larger compared to $y(x)$ but overall they both have not too small error.

  • 1
    $\begingroup$ Is your data just the solution you get from your finite difference method? Then your finite difference approximation of y' would be your best approximation based on the data. Then the cubic spline approach you use is likely the best you can do as you only have second order accuracy in your solution method. I'm not familiar with NDSolve's exact algorithms, but I suspect is of higher order of accuracy, and likely uses more nodes than your non-adaptive method. $\endgroup$
    – user9794
    Dec 10, 2022 at 19:02
  • $\begingroup$ Have you considered reducing this to a system of two first-order equations (one of which is an equation for $y'$)? $\endgroup$ Dec 10, 2022 at 22:35
  • $\begingroup$ @user9794 Yes, the data I got is from the finite difference method. From what you said, I presume that the error I have for y' based on the interpolation function of y is quite normal. I guess what I'm thinking is if there are ways to obtain y' with less error. I'm quite new to numerical analysis so I'm not so familiar with other methods or any backdoor techniques. $\endgroup$
    – mathemania
    Dec 11, 2022 at 7:34
  • $\begingroup$ @BrianBorchers Not as of the moment, I haven't thought of doing that since I'm just doing all calculations based on how numerical analysis textbooks are doing it, and I haven't seen any finite difference method where the ODE is split into two first order equations. $\endgroup$
    – mathemania
    Dec 11, 2022 at 7:36

2 Answers 2


In my version of Mathematica it appears that the NDSolve solution is returning a cubic Hermite spline, though I have no idea what order it is using internally to solve. One thing which is immediately obvious from plotting the residual of your NDSolveValue solution is the number of points being used is much higher than 10, so even if it was 2nd order accurate like your method it shouldn't be surprising that it is more accurate overall.

Another potential source of error is that I don't know what type of boundary conditions Mathematica is using for the spline interpolation. I get yfunc''[0]=-0.500231, which certainly doesn't satisfy the original ODE since you're supposed to have $y''(0) = 0$. I don't know if there's any way to tell Mathematica to change interpolation boundary conditions.

You can get an overall picture of how good your (interpolated) solution is by plotting the residual:

Plot[{yfunc''[x] + 3 yfunc'[x] yfunc[x]}, {x, 0, 2}, Frame -> True, PlotRange -> All]

enter image description here

Other than doing interpolation and taking the derivative of the interpolating function, you can do something similar to what your Mathematica implementation is doing and split your ODE into a coupled system of first order ODE's and simultaneously solve for $y$ and $y'$:

$$ v' + 3 y v = 0\\ y' = v\\ y(0) = 0\\ y(2) = 1 $$ The downside of this is I think you will also have to use a shooting method since $v(0)$ is unknown. It's also not necessarily more accurate since if you use second order finite difference approximations you'll still get at best a second order solution (might even end up with $y'$ being first order while $y$ is second order, I haven't done the error analysis of this method).


I do find it strange that you are decomposing your ODE into a coupled system of first order ODEs then doing a shooting method. NDSolveValue can solve boundary value problems:

sols2 = NDSolveValue[{y''[x] == (-3 y[x] y'[x]), y[0] == 0, y[2] == 1}, y, {x, 0, 2}, AccuracyGoal -> 20, PrecisionGoal -> 20, WorkingPrecision -> 35, MaxSteps -> Infinity];
  • $\begingroup$ Nice analysis, so it means that if I'm interested up to y' only (to be used in some other things) then my data for the derivative is kinda ok? I mean, y'' doesn't seem to be correct to a large tolerance but the error of y' is still kinda acceptable, I guess. $\endgroup$
    – mathemania
    Dec 12, 2022 at 6:07
  • $\begingroup$ I suppose it depends on what you consider acceptable, though there's not anything fundamentally wrong with your approach other than not satisfying $y''(0)$, which may or may not matter. You can of course always increase the number of points you use for the finite difference method and it should become more accurate. $\endgroup$ Dec 13, 2022 at 3:44

Here is another, more analytical view on your problem: You can also rewrite your equation as $$ y^{\prime\prime} + \frac32 (y^2)^\prime = 0 $$ and thus reduce one derivative term: $$ y^{\prime} + \frac32 y^2 = C \tag{1}\label{1} $$

Using the left boundary conditions, one can calculate the constant $C$: $$ C = y^{\prime}|_{x=0} $$

Basically, Eq. (\ref{1}) gives you a recipe for calculating the derivative in terms of $y$. However, for this the constant $C$ is required, and the (one-sided) derivative condition in the previous equation will likely not not be too accurate.

Therefore, using separation of variables one can solve the equation analytically:

$$ \frac{dy}{C - \frac32 y^2} = dt\\ \rightarrow y = \sqrt \frac{2}{3C} \text{arctanh}\left(\sqrt{ \frac 3{2C}} \ x\right) + C_2 $$

From the left boundary condition, $y(0) = 0$, one directly obtains $C_2 = 0$.

Calculating the constant $C$ from the second boundary condition is more difficult, and I'm not aware of an analytical solution. However, you can solve for $C$ via Newton iteration, which then let's you tackle the derivative via Eq. (\ref{1}).

But frankly, then you also have the sought solution $y$ and probably won't need the finite-difference solution at all.


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.