2
$\begingroup$

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]

Mathematica graphics Mathematica graphics Mathematica graphics

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

$\endgroup$
  • $\begingroup$ @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. $\endgroup$ – Hugh Mar 2 '18 at 9:25
  • $\begingroup$ Yes, the accuracy of numerical sulutions of any finite differences or finite element solver should reduce for higher derivatives. $\endgroup$ – Henrik Schumacher Mar 2 '18 at 9:27
1
$\begingroup$

A possible solution is to choose a higher difforder. The following is obtained with points = 100; difforder = 8:

Mathematica graphics

The troublesome part is, FindRoot becomes slow when points get larger. (That's why I choose points = 100 here. )

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Hugh Mar 2 '18 at 11:59
  • $\begingroup$ @hugh Slightly different. Your approach is amount to raising the difference order only for the high order derivatives. This might make the calculation faster because the resulting system is less intricated. If you want to have a try, you may find this post interesting. $\endgroup$ – xzczd Mar 2 '18 at 12:17
  • $\begingroup$ Interesting post. I have just managed to use all first order equations but getting more points is making everything very slow. Another point. If I use first order equations which equations should I drop? Equations at the start or equations at the end of the grid? Does it make much difference? Thanks again for your help. $\endgroup$ – Hugh Mar 2 '18 at 15:53
  • $\begingroup$ @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)$. $\endgroup$ – xzczd Mar 3 '18 at 3:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.