# Curve Fitting ODE

Parameters

n = 0.6;
Re1 = 1;
Re2 = 1;
m = 2;


The system of the differential equations

soln1 = NDSolve[{Re1^2 u[r]/r == p'[r],
u''[r] + u'[r]/ r -  u[r]/r^2 - 2 n/Re2 w'[r] == 0,
u'[r] + u[r]/ r + (2 - n/(m^2 Re2)) (w''[r] + w'[r]/ r ) +
4 Re2 w[r] == 0, u[0.01] == 0, u[1] == 1, w[0.01] == 0,
w[1] == 1, p[0.01] == 0}, {u, w, p}, {r, 0.01, 1}];


Plotting

uF = u[r] /. soln1;
wF = w[r] /. soln1;
pF = p[r] /. soln1;
Plot[{uF}, {r, 0.01, 1}, LabelStyle -> 20, Frame -> True,
FrameLabel -> {{" u", "  "}, {" r", " "}}, PlotStyle -> Blue,
Axes -> False]
Plot[{wF}, {r, 0.01, 1}, LabelStyle -> 20, Frame -> True,
FrameLabel -> {{" w", "  "}, {" r", " "}}, PlotStyle -> Blue,
Axes -> False]
Plot[{pF}, {r, 0.01, 1}, LabelStyle -> 20, Frame -> True,
FrameLabel -> {{" p", "  "}, {" r", " "}}, PlotStyle -> Blue,
Axes -> False]


I want to do a curve fitting to find a expression of polynomial x and x^2 for uF, wF and pF

Fit[uF, {1, r, r^2}, r]

Plot[%, {r, 0.01, 1}]


I need it to match with uF, wF and pF as closer as possible

I think the range of parameters is the essential thing to reach this goal

what do you think?

• "As close as possible": With respect to what norm? Jun 16 '17 at 17:37
• close to numerical solution that I gain from NDsolving set of ODE (uF, wF and pF) Jun 16 '17 at 22:32
• You usually measure how close functions with some sort of norm. Common norms are the infinity norm, 1-norm, and 2-norm (which is similar to the root-mean-square error). The 1-norm and 2-norm are defined in terms of an integral. The infinity norm is equal to the maximum error over the interval. One also can measure the relative error instead of the absolute error. (Strictly speaking one cannot use relative error if the function is zero on the interval, as in your case.) Jun 16 '17 at 23:35

## The "FunctionApproximations" package

There is a "FunctionApproximations" package that comes with Mathematica that seems designed for this sort of problem. The functions MiniMaxApproximation and GeneralMiniMaxApproximation can find a minimax approximant over a finite interval with respect to relative or absolute error (or in general a weight function). If the function (or weight function) being approximated crosses zero, as in the OP, then the approximant with respect to relative error is impossible to achieve. If the weight function never changes sign, one can take an approach similar to that of AccuracyGoal and ScaledVectorNorm, and add a little bias or offset to the function or weight function before carrying out the approximation. That said, the minimax algorithm is finicky, and one often has to adjust the option settings of the options Brake, Bias, and MaxIterations.

## The OP's setup

There are a few issues in the OP's problem. If you are going to be approximating a solution obtained from NDSolve, I recommend using the option InterpolatingOrder -> All, unless it is prohibitively expensive. This minimizes the interpolation error between the steps. Further, approximation methods are usually sensitive to singularities (in the function or its derivatives). An InterpolatingFunction generally has a singularity at each step. So I chose the integration method "Extrapolation" to reduce the number of steps. Finally, one sometimes needs a higher working precision to conquer numerical instability. I'll present only solutions that do not need it, but for some degrees, one needs higher than machine precision.

ClearAll[u, r, p, q, s, t, w, x, y, z, a, b, c, m, n];
n = 6/10;
Re1 = 1;
Re2 = 1;
m = 2;

soln1 = NDSolve[{Re1^2 u[r]/r == p'[r],
u''[r] + u'[r]/r - u[r]/r^2 - 2 n/Re2 w'[r] == 0,
u'[r] + u[r]/r + (2 - n/(m^2 Re2)) (w''[r] + w'[r]/r) +
4 Re2 w[r] == 0, u[0.01] == 0, u[1] == 1, w[0.01] == 0,
w[1] == 1, p[0.01] == 0}, {u, w, p}, {r, 0.01, 1},
Method -> "Extrapolation", InterpolationOrder -> All(*, WorkingPrecision -> 32*)];

uF = u /. First@soln1;
wF = w /. First@soln1;
pF = p /. First@soln1;


## Solutions

I'll look at rational approximations. The function wF is almost vertical at r == 0.01 and quickly turns over. It's difficult to approximate with a low-degree polynomial (the Chebyshev series of degree 16 has a similar accuracy to the degree-(3,3) minimax rational approximation shown below.)

Needs["FunctionApproximations"];


### Minimizing relative error

The degree 2 over degree 3 approximant wApprox does a bit better than a least squares fit as in MikeY's approach. It is necessary to raise MaxInterations to get convergence. Here we add the bias acc (acc for "accuracy," as in AccuracyGoal) to wF[r] and subtract it from the approximant returned by MiniMaxApproximation. One can see the characteristic equioscillation of the error of the minimax approximation in the relative error plot.

(* Minimax (wF[r] - approx[r])/(acc + wF[r]) *)
acc = 1*^-8; (* cf. AccuracyGoal -> 8 *)
approx = MiniMaxApproximation[acc + wF[r], {r, {0.01, 1}, 2, 3},
MaxIterations -> 200];
{wApprox, relerror} =             (* see code dump at end for an explanation of Replace *)
Replace[approx, {_List, {ap_, e_}} | {_List, ap_, e_} :> {ap, Max@Abs@e}]
(*
{(-0.69499 + 68.6614 r + 83.7666 r^2) /
(1 + 67.268 r - 63.949 r^2 + 143.785 r^3),          (* approximant of wF *)
0.0245034}                                            (* error *)
*)

plotMMA[wF[r], approx, acc, {r, 0.01, 1}]


The following uses GeneralMiniMaxApproximation to implement the bias with a weight function acc + wF[t] for the degree-(3,3) approximant.

(* Minimax (wF[r] - approx[r])/(acc + wF[r]) *)
acc = 1*^-8;
approx =
GeneralMiniMaxApproximation[{t, wF[t], acc + wF[t]}, {t, {0.01, 1}, 3, 3}, r,
MaxIterations -> 200];
{wApprox, abserror} =
Replace[approx, {_List, {ap_, e_}} | {_List, ap_, e_} :> {ap, Max@Abs@e}]
(*
{(-0.849961 + 73.5664 r + 1151.73 r^2 - 875.442 r^3) /
(1 + 119.214 r + 461.739 r^2 - 231.094 r^3),        (* approximant of wF *)
0.00529828}                                           (* error *)
*)

plotMMA[wF[r], approx, {0, acc}, {r, 0.01, 1}]


### Minimizing the absolute error

To minimize the absolute error, use the weight function 1. (No need to use an offset acc.) Here, to get convergence, we needed to tweak Brake and Bias. Increasing Brake (the default is {5, 5}) slows down the steps. Sometimes the algorithm steps out of the feasible region (which is not guaranteed to exist for all cases), and Brake can help rein it in. Bias skews the the initial interpolation abscissas. From the plot of wF[r], we can see we might need denser sampling on the left near r == 0.1, where the graph hits the axis like a square root. Setting Bias to a negative value skews the sampling to the left. (Setting it less than -0.5 seems to make the sampling start outside r == 0.01, which fails because wF is an InterpolatingFunction with domain {0.01, 1.}).

(* Minimax (wF[r] - approx[r]) *)
approxAbs =
GeneralMiniMaxApproximation[{t, wF[t], 1}, {t, {0.01, 1}, 3, 3}, r,
Brake -> {100, 30}, Bias -> -0.5, MaxIterations -> 200];
{wApproxAbs, abserror} =
Replace[approxAbs, {_List, {ap_, e_}} | {_List, ap_, e_} :> {ap, Max@Abs@e}]
(*
{(-0.664799 + 62.2947 r + 464.679 r^2 - 350.072 r^3) /
(1 + 74.7161 r + 154.849 r^2 - 53.9008 r^3),           (* approximant of wF *)
0.00242014}                                             (* error *)
*)

plotMMA[wF[r], approxAbs, 0, {r, 0.01, 1}]


### Degree-(3,2) is a difficult case

It seems difficult to compute the degree-(3,2) approximant to wF[r]. Indeed, with an odd-degree numerator and a degree-2 denominator, one gets a spike in the error around 0.01 < r < 0.085. If one splits the interval at r == 0.15, one can get excellent approximants on each interval. (I was trying to compare a degree-2 denominator solution to MikeY's. It's similar in error. But really, MiniMaxApproximation failed to find one of the extrema, the greatest one at that; and it should have reported an error, MiniMaxApproximation::extalt, too many extrema.)

(* Minimax (wF[r] - approx[r])/(acc + wF[r]) *)
acc = 1*^-8; (* cf. AccuracyGoal -> 8 *)
approx = MiniMaxApproximation[acc + wF[r], {r, {0.01, 1}, 3, 2},
MaxIterations -> 200, Brake -> {30, 30}];
{wApprox, relerror} =
Replace[approx, {_List, {ap_, e_}} | {_List, ap_, e_} :> {ap, Max@Abs@e}]
(*
{(-0.471551 + 45.8908 r + 127.465 r^2 - 103.279 r^3) /
(1 + 40.4979 r + 28.1309 r^2),                      (* approximant of wF *)
0.000332984}                                          (* error *)
*)
plotMMA[wF[r], approx, acc, {r, 0.01, 1}]


### A high-precision solution

The solution by NDSolve has a precision goal of about 8, so pursuing a better approximant than this does not make a lot of sense.

(* Minimax (wF[r] - approx[r])/(acc + wF[r]) *)
acc = 1*^-8;
approx =
GeneralMiniMaxApproximation[{t, wF[t], acc + wF[t]}, {t, {0.01, 1}, 10, 6}, r,
Brake -> {100, 30}, Bias -> -0.5, MaxIterations -> 200];
{wApprox, abserror} =
Replace[approx, {_List, {ap_, e_}} | {_List, ap_, e_} :> {ap, Max@Abs@e}]
(*
{(-1.58409 - 521.421 r + 19534.7 r^2 + 4.03048*10^6 r^3 + 7.84636*10^7 r^4 +
2.97461*10^8 r^5 - 2.98576*10^7 r^6 - 2.51189*10^8 r^7 + 4.20181*10^7 r^8 +
4.5839*10^7 r^9 - 1.27469*10^7 r^10) /
(1. + 988.287 r + 151504. r^2 + 5.84738*10^6 r^3 + 6.19884*10^7 r^4 +
1.43813*10^8 r^5 - 3.7763*10^7 r^6),                (* approximant of wF *)
1.68198*10^-8}                                         (* error *)
*)

plotMMA[wF[r], approx, {0, acc}, {r, 0.01, 1}]


## plotMMA code

There are two forms of return value of MiniMaxApproxition, depending on whether it was successful or almost successful (MiniMaxApproximation::extalt trouble with the extrema or MaxIterations):

{List of abscissae of error extrema, {approximant, error}}
{List of abscissae of error extrema, approximant, List of error extrema}


Thus Replace[approx, {_List, {ap_, e_}} | {_List, ap_, e_} :> ap] extracts the approximant function from approx. Above in the examples, the following was used, which extracts the approximant and maximum error:

Replace[approx, {_List, {ap_, e_}} | {_List, ap_, e_} :> {ap, Max@Abs@e}]


Plot routine:

plotMMA[fn_, approx_, bias_, {x_, a_, b_}] := plotMMA[fn, approx, {bias, bias}, {x, a, b}];
plotMMA[fn_, approx_, {bias_, acc_}, {x_, a_, b_}] :=
Module[{degrees, wApprox},
wApprox = Replace[approx, {_List, {ap_, e_}} | {_List, ap_, e_} :> ap];
degrees = Exponent[Through[{Numerator, Denominator}[wApprox]], r];
wApprox = wApprox - bias;
GraphicsRow[
Plot[#1, {x, a, b},
PlotRange -> All, PlotLabel -> #2,
GridLines -> {approx[[1]], None}] &,
{
{1/(acc + wF[r]) (wF[r] - (wApprox)),
wF[r] - (wApprox),
{wF[r], wApprox}},
{"Relative error", "Error", "Function, approximation"}
}
],
PlotLabel -> Row[{"Degrees ", degrees}]
]
];

• As an OBTW, if you take my solution for wF and turn it into a rational approximation, it ends up being Degree(5,2). Jun 19 '17 at 16:42
• @MikeY Yes. OTOH, your least-squares fit has 6 degrees of freedom, which is more like degree (3,2) or (2,3). It seemed more appropriate to make that comparison. Jun 19 '17 at 19:29

First, the Fit function needs data to fit to, so you need to take your uF and generate some data. This is coarse just to show how.

dat = Table[{r, uF // First}, {r, 0.01, 1, .05}]

(*   {{0.01, 0.}, {0.06, 0.0335532}, {0.11, 0.0788975}, {0.16,
0.129989}, {0.21, 0.184636}, {0.26, 0.241593}, {0.31,
0.299998}, {0.36, 0.359177}, {0.41, 0.418568}, {0.46,
0.47768}, {0.51, 0.536074}, {0.56, 0.593351}, {0.61,
0.649142}, {0.66, 0.703109}, {0.71, 0.754937}, {0.76,
0.804337}, {0.81, 0.851044}, {0.86, 0.894815}, {0.91,
0.935435}, {0.96, 0.972707}} *)


Then you can fit the data. Eyeballing the curve of uF, you can see that it has a change in the 2nd derivative, so you need at least a 3rd order function.

cf = Fit[dat, {1, r, r^2, r^3}, r]

(*  -0.0152527 + 0.810384 r + 0.860931 r^2 - 0.663355 r^3 *)

Plot[{cf, uF}, {r, 0.01, 1}]


Edit: for the wF function.

To approximate the wF curve, in looking at it, you need a little something that dissipates with r. A typical function is 1/r or 1/r^2. Do this...

dat = Table[{r, wF // First}, {r, 0.01, 1, .01}];

cf = Fit[dat, {1, r, 1/r, 1/r^2, r^2, r^3}, r];

Plot[{wF // First, cf}, {r, 0.01, 1}, PlotLegends -> {"original", "fitted"}]
`

• thank u so much MikeY , it work for uF and pF but for wF didn't match exactly and I do not want to increase the order of function Jun 16 '17 at 15:38
• any suggestion @ MikeY Jun 16 '17 at 15:43
• See my edit to the answer, you have to add a couple of terms to get a great fit. Jun 16 '17 at 16:07