# Using NIntegrate in Finite Difference Derivative method

I'm trying to solve a second order differential equation using the code given by @xzczd here which is based on this.

What this code makes is transform differential equations in algebraic ones by means of the Finite Difference method.

When I solve the equation

    R*X''[R] - X'[R] + R^3 == 0;


I have no problems at all. However when I change the equation to

R*X''[R] - X'[R] + R^3*nColdr[R] == 0;


being that

nColdr[R_?NumberQ]=NIntegrate[g[x],{x,0,R}],


where g is a complicated function of x, I get the output:

FindRoot::nlnum: The function value {0. +2.09232*10^17 (-1.37456*10^-14+2.7573*10^-14 nColdr$6323491[{0.,0.20202,0.40404,0.606061,0.808081,1.0101,1.21212,1.41414,1.61616,1.81818,2.0202,2.22222,2.42424,<<26>>,7.87879,8.08081,8.28283,8.48485,8.68687,8.88889,9.09091,9.29293,9.49495,9.69697,9.89899,<<50>>}]),<<49>>,<<150>>} is not a list of numbers with dimensions {200} at {X$6323491[0],X$6323491[20/99],X$6323491[40/99],X$6323491[20/33],X$6323491[80/99],X$6323491[100/99],X$6323491[40/33],X$6323491[140/99],X$6323491[160/99],X$6323491[20/11],X$6323491[200/99],X$6323491[20/9],X$6323491[80/33],X$6323491[260/99],X$6323491[280/99],X$6323491[100/33],X$6323491[320/99],<<17>>,X$6323491[680/99],X$6323491[700/99],X$6323491[80/11],X$6323491[740/99],X$6323491[760/99],X$6323491[260/33],X$6323491[800/99],X$6323491[820/99],X$6323491[280/33],X$6323491[860/99],X$6323491[80/9],X$6323491[100/11],X$6323491[920/99],X$6323491[940/99],X$6323491[320/33],X$6323491[980/99],<<150>>} = {1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,<<150>>}. >>


I'm pretty sure that the problem is because of the NIntegrate... I already try to do Hold and then ReleaseHold in different parts of the @xzczd code but it's not working...

It's because pdetoode/pdetoae has made use of Listable attribute, which almost all the arithmetic functions own. So you just need to make your nColdr Listable. A simple example:

nColdr[x_?NumberQ] := NIntegrate[y^2, {y, 0, x}]

SetAttributes[nColdr, Listable]

eq = R*X''[R] - X'[R] + R^3 nColdr[R] == 0;

bc = {X[0] == 0, X[1] == 1};

domain = {0, 1};
points = 25;
grid = Array[# &, points, domain];
difforder = 4;
ptoafunc = pdetoae[X[R], grid, difforder];

ae = Delete[ptoafunc[eq], {{1}, {-1}}];

data = FindRoot[{ae, bc}, {X@#, 0} & /@ grid][[All, -1]];

pic = ListLinePlot[data, DataRange -> domain];

eqref = R*X''[R] - X'[R] + R^3 Integrate[y^2, {y, 0, R}] == 0;
sol = NDSolveValue[{eqref, bc /. X[0] -> X[10^-6]}, X, {R, 10^-6, 1}];
ListLinePlot[sol, PlotStyle -> {Red, Dashed, Thick}]~Show~pic


• Thank you! Let me just ask the following: I thought that the variable 'difforder' was supposed to indicate the order of the differential equation, but now I see I was completely wrong since you set it to 4. What is then the meaning of it? – AJHC Jun 6 '18 at 15:59
• @AJHC It means "difference order" i.e. the order of difference formula being chosen. (Perhaps I should not use this abbreviation any more from now on… ) – xzczd Jun 6 '18 at 16:00
• It's okay I think... the problem is mine because I couldn't understand the pdetoode code. Anyway, as higher the difference order the higher the precision right? – AJHC Jun 6 '18 at 16:35
• @AJHC Nope, this is a rather complicated topic actually. Some discussion can be found in tutorial/NDSolvePDE of the document. A rule of thumb is, 4 or 2 will work for most cases, and "Try it out" is the only way to figure out the best suited difference order in practice. – xzczd Jun 6 '18 at 16:43
• @ajhc It depends on the position of the other boundary condition, usually the difference equations that are the closest ones to the b.c. are to be removed. If it's e.g. X[0]==… then ae = Delete[ptoafunc[eq], {{1}, {2}}];, if it's e.g. X[1]==… then ae = Delete[ptoafunc[eq], {{1}, {-1}}];. Your next question is probably "why I should remove difference equation in this way?" Well, to be honest, I don't know, this is somewhat a empirical rule to me. – xzczd Jul 10 '18 at 13:10