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


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...


1 Answer 1


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

Mathematica graphics

  • $\begingroup$ 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? $\endgroup$
    – AJHC
    Commented Jun 6, 2018 at 15:59
  • $\begingroup$ @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… ) $\endgroup$
    – xzczd
    Commented Jun 6, 2018 at 16:00
  • $\begingroup$ 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? $\endgroup$
    – AJHC
    Commented Jun 6, 2018 at 16:35
  • $\begingroup$ @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. $\endgroup$
    – xzczd
    Commented Jun 6, 2018 at 16:43
  • 1
    $\begingroup$ @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. $\endgroup$
    – xzczd
    Commented Jul 10, 2018 at 13:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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