2

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

xzczd
  • 65,995
  • 9
  • 163
  • 468
AJHC
  • 369
  • 1
  • 8

1 Answers1

3

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

xzczd
  • 65,995
  • 9
  • 163
  • 468
  • 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 06 '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 06 '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 06 '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 06 '18 at 16:43
  • Maybe the tutorial you are refering to is tutorial/NDSolveMethodOfLines instead? – AJHC Jun 07 '18 at 08:37
  • @AJHC Oh, I forgot its name has changed after v10. Anyway, I'm referring to the tutorial that involves introduction for NDSolve`FiniteDifferenceDerivative. – xzczd Jun 07 '18 at 10:09
  • I can see in the tutorial that for each difference order there are, at least, three different formulas for $f'$ and $f''$. Is it possible to know which one Mathematica is using? – AJHC Jun 07 '18 at 20:59
  • @AJHC NDSolve`FiniteDifferenceDerivative uses one-sided difference formula at and near the boundary, and central difference formula in the middle, or it won't be able to generate enough equations for the whole domain. You can use pdetoae to discretize a single X'[R] and observe the output. If you're not familiar with one-sided formula, start from page 6 of this book. – xzczd Jun 08 '18 at 06:10
  • I'm now trying to use pdetoae with the BCs X'[0]==X[0] and I also tried X'[0]==0 but in both cases I get the error The function value etc etc is not a list of numbers with dimensions...

    Isn't it possible to use pdetoae with Neuman-type BCs?

     – AJHC Jul 01 '18 at 14:40
  • @ajhc Yes, just try it out. – xzczd Jul 01 '18 at 14:50
  • For example, using exactly the same code of your answer, but using X'[0]==0 instead of X[0]==0 shouts the error The function value (...) is not a list of numbers with dimensions {25} at (...) – AJHC Jul 01 '18 at 15:01
  • @AJHC That's because the b.c.s in your question are Dirichlet b.c.s so there's no need to use pdetoae on them and I did not. Now you just need to use pdetoae on your b.c., try it out. – xzczd Jul 01 '18 at 15:04
  • you're right. I changed the condition X'[0]==X[0] from a BC to an equation which made me stay with only one BC. Then I had 10+1 algebraic equations, where the first one from the group of 10 and the new one were incompatible, and so I deleted the first one from the group of 10. However, I still had to delete one more equation from the group of 10 and I chose it randomly... Is that ok? – AJHC Jul 10 '18 at 12:32
  • @ajhc …Discretization of X'[0]==X[0] doesn't need to be that troublesome, try ptoafunc[X'[0]==X[0]]. – xzczd Jul 10 '18 at 12:39
  • Yes, that gives the same thing with less effort ;) my question is: when deleting 2 of the equations, one of them I delete because it's incompaible with the equation that results from ptoafunc[X'[0]==X[0]]. What about the other one? I can chose it randomly? – AJHC Jul 10 '18 at 12:50
  • 1
    @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