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.)
First load the package:
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[
MapThread[
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}]
]
];
