2
$\begingroup$

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

Improve this question
$\endgroup$
3
$\begingroup$

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

Improve this answer
$\endgroup$
16
  • $\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 Jun 6 '18 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 Jun 6 '18 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 Jun 6 '18 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 Jun 6 '18 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 Jul 10 '18 at 13: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.