# Oscillations on solution of finite difference equation

In a previous post on the solution of an ODE with a boundary conditon at infinty I had some excelent help from xzczd and am now returning with a further problem along the same lines. I have used the code pdetoae which may be found in this link.

The problem is that the solution has oscillations. These oscillations are particularly clear on the first and second derivatives of the solution which I need.

The ODE, initial and final conditions and the variables are:

    inf = 15;
nn = 2; eqns = {(-α +
2*α*Derivative[1][a[0]][y]^2 + α*
Derivative[1][a[1]][y]^2 + α*
Derivative[1][a[2]][y]^2 + α*
Derivative[1][b[1]][y]^2 + α*Derivative[1][b[2]][y]^2 -
2*α*a[0][y]*Derivative[2][a[0]][y] - α*a[1][y]*
Derivative[2][a[1]][y] - α*a[2][y]*
Derivative[2][a[2]][y] - α*b[1][y]*
Derivative[2][b[1]][y] - α*b[2][y]*
Derivative[2][b[2]][y] - 2*Derivative[3][a[0]][y])/2 ==
0, (4*α*Derivative[1][a[0]][y]*Derivative[1][a[1]][y] +
2*α*Derivative[1][a[1]][y]*Derivative[1][a[2]][y] +
2*Derivative[1][b[1]][y] +
2*α*Derivative[1][b[1]][y]*Derivative[1][b[2]][y] -
2*α*a[1][y]*Derivative[2][a[0]][y] -
2*α*a[0][y]*Derivative[2][a[1]][y] - α*a[2][y]*
Derivative[2][a[1]][y] - α*a[1][y]*
Derivative[2][a[2]][y] - α*b[2][y]*
Derivative[2][b[1]][y] - α*b[1][y]*
Derivative[2][b[2]][y] - 2*Derivative[3][a[1]][y])/2 ==
0, (2 - 2*Derivative[1][a[1]][y] +
4*α*Derivative[1][a[0]][y]*Derivative[1][b[1]][y] -
2*α*Derivative[1][a[2]][y]*Derivative[1][b[1]][y] +
2*α*Derivative[1][a[1]][y]*Derivative[1][b[2]][y] -
2*α*b[1][y]*Derivative[2][a[0]][y] - α*b[2][y]*
Derivative[2][a[1]][y] + α*b[1][y]*
Derivative[2][a[2]][y] -
2*α*a[0][y]*Derivative[2][b[1]][y] + α*a[2][y]*
Derivative[2][b[1]][y] - α*a[1][y]*
Derivative[2][b[2]][y] - 2*Derivative[3][b[1]][y])/2 ==
0, (-α + α*Derivative[1][a[1]][y]^2 +
4*α*Derivative[1][a[0]][y]*
Derivative[1][a[2]][y] - α*Derivative[1][b[1]][y]^2 +
4*Derivative[1][b[2]][y] -
2*α*a[2][y]*Derivative[2][a[0]][y] - α*a[1][y]*
Derivative[2][a[1]][y] -
2*α*a[0][y]*Derivative[2][a[2]][y] + α*b[1][y]*
Derivative[2][b[1]][y] - 2*Derivative[3][a[2]][y])/2 ==
0, (-4*Derivative[1][a[2]][y] +
2*α*Derivative[1][a[1]][y]*Derivative[1][b[1]][y] +
4*α*Derivative[1][a[0]][y]*Derivative[1][b[2]][y] -
2*α*b[2][y]*Derivative[2][a[0]][y] - α*b[1][y]*
Derivative[2][a[1]][y] - α*a[1][y]*
Derivative[2][b[1]][y] -
2*α*a[0][y]*Derivative[2][b[2]][y] -
2*Derivative[3][b[2]][y])/2 == 0};
ic0 = {a[0][0] == 0, Derivative[1][a[0]][0] == 0, a[1][0] == 0,
Derivative[1][a[1]][0] == 0, b[1][0] == 0,
Derivative[1][b[1]][0] == 0, a[2][0] == 0,
Derivative[1][a[2]][0] == 0, b[2][0] == 0,
Derivative[1][b[2]][0] == 0};
bc = {Derivative[2][a[0]][15] == 0, Derivative[1][a[1]][15] == 1,
Derivative[1][b[1]][15] == 0, Derivative[1][a[2]][15] == 0,
Derivative[1][b[2]][15] == 0};
var = {a[0], a[1], b[1], a[2], b[2]};


The code to make the grid, use pdetoae and solve is:

points = 200;
grid = Array[# &, points, {0, inf}];
difforder = 4;
ptoafunc = pdetoae[var[y], grid, difforder];
del = #[[3 ;; -2]] &;
ae = del /@ ptoafunc@eqns;
aebc = ptoafunc@{ic0, bc};
sol = Partition[
With[{guess = 1},
FindRoot[{ae, aebc} /. {α -> 0.5},
Flatten[Table[{var[[i]][y], guess}, {i, 2 nn + 1}, {y, grid}],
1],
MaxIterations -> 500]][[All, -1]], points];


I now look at the solutions and their first and second differences:

diffs = Differences[#] & /@ sol;
diffs2 = Differences[#, 2] & /@ sol;
sols = Transpose[{grid, #}] & /@ sol;
dsols = Transpose[{Most[grid], #}] & /@ diffs;
dsols2 = Transpose[{Drop[grid, -2], #}] & /@ diffs2;
ListLinePlot[sols, PlotRange -> All]
ListLinePlot[dsols, PlotRange -> All]
ListLinePlot[dsols2, PlotRange -> All]


As you can see the solution looks good but the first and second differences show an oscillation that is associated with the grid. Is this inevitable with a finite difference method?

I can smooth the solution using various methods but I am not sure if this is valid. How do I get rid of the oscillations?

Thanks

• @xzczd Sorry to bother you again with further questions on your excellent method. However, if you could suggest a reason for these oscillations and a method for dealing with them that would help. – Hugh Mar 2 '18 at 9:25
• Yes, the accuracy of numerical sulutions of any finite differences or finite element solver should reduce for higher derivatives. – Henrik Schumacher Mar 2 '18 at 9:27

A possible solution is to choose a higher difforder. The following is obtained with points = 100; difforder = 8:
The troublesome part is, FindRoot becomes slow when points get larger. (That's why I choose points = 100 here. )
• Thanks for the fast response. Would it help if I reduced the order (here three) of the equations by introducing three first order equations? Does it come to the same thing in the end as a problem for FindRoot? – Hugh Mar 2 '18 at 11:59
• @hugh Well, I should admit this is something I didn't explore much. In principle you can just use #[[3;;-2]]& on the equations generated from the orignal equations, but it's probably not the best way because one-sided formula is used for calculating derivative values on the boundary in this case. I think the best way is the following: use e.g. $y_1(0)=1$ (here $y_1$ represents $y'$) to replace the left boundary of auxiliary ODE $y_1'(x)=y_2(x)$. – xzczd Mar 3 '18 at 3:10