I am solving PDEs with NDSolve and NIntegrate, but I do not know how to pass arguments correctly. My orignial code is very complicated, so I simplfied it to clarify the problems I met.

1) Firstly, I solved the simplified code without passing the arguments:

ClearAll[x, y,r, t];
ss = NDSolve[{
     x'[t] == (NIntegrate[x[t]*r*r, {r, 0, 5}] + y[t] + t)*y[t], 
     y'[t] == -x[t], x[0] == 1, y[0] == 1}, {x, y}, {t, 10}];
Plot[Evaluate[{x[t], y[t]} /. ss], {t, 0, 10}]

2) Then I used the skil of passing arguments to solve the problems, and the code becomes:

ClearAll[x, y, z, zz, t];
u[x_, r_] := x*r;
v[u_?NumericQ] := NIntegrate[u*r, {r, 0, 5}];
z[v_, y_, t_] := v + y + t;

s = NDSolve[{
    x'[t] == z[t, t, t]*y[t], 
    y'[t] == -x[t], x[0] == 1, y[0] == 1}, {x, y}, {t, 10}];
Plot[Evaluate[{x[t], y[t]} /. s], {t, 0, 10}]

But I got different and wrong results. Why?

3) I also tried to pass the arguments in another way, but I just got the error messages. Why?

ClearAll[x, y, u, v, z, zz, t, f];
u[x_, r_] := x*r;
v[x_?NumberQ] := NIntegrate[u*r, {r, 0, 5}];
z[v_, y_, t_] := v + y + t;

s = NDSolve[{
    x'[t] == z[v[u[x[t]]], y[t], t]*y[t], 
    y'[t] == -x[t], x[0] == 1, y[0] == 1}, {x, y}, {t, 10}];

4) I also tried another way to pass the arguments, and the result is correct and same to 1). But I prefer not to use it due to my complicated original code.

ClearAll[x, y, u, v, z, zz, t, f];
u[x_, r_] := x*r;
v[x_?NumberQ] := NIntegrate[u*r, {r, 0, 5}];
f = NIntegrate[u[x[t], r]*r, {r, 0, 5}];
z[v_, y_, t_] := v + y + t;

s = NDSolve[{x'[t] == z[f, y[t], t]*y[t], y'[t] == -x[t], x[0] == 1, 
y[0] == 1}, {x, y}, {t, 10}];

Plot[Evaluate[{x[t], y[t]} /. s], {t, 0, 10}]

How to pass arguments in complicated codes correctly?

The original equations look like this. They are integration-differential equations.

The only independent variabls is z, the aim is to solve $\nu$[z], qr[z] and qi[z] from z=0 to z. The parameters are: A[z], B[z], D[z], d[z], and ne[$\rho$,z], $\gamma$[$\rho$,z], $\epsilon$r[$\rho$,z] and $\epsilon$i[$\rho$,z]. The integration of $\rho$ ranges from d to infinite.


enter image description here

  • * I wrote the code with error messages: NDSolve::ndnum: Encountered non-numerical value for a derivative at z == 0.1.

From equations to code, different names are used: $\rho$ -->r, D-->D0, d-->dmin. $\gamma$ -->$\gamma_0$

ClearAll[ν, qr, qi, r, γ, ne, ϵr, ϵi, dmin, A, B, D0, aL, Zni, νei];
aL = 10;
Zni = 0.001;
νei = 0.001;
γ= (1 + aL^2*ν^2*Exp[-2*qr*r^2])^0.5;
ne= Zni - 4*qr*γ*(1 -γ^(-2))*(1 - qr*r^2*(1 + γ^(-2)));
ϵr= 1 - ne/γ;
ϵi= 4*Pi*νei*ne;
dmin= Re[qr^(-0.5)*(-0.5*Log[Zni/(8*aL^2*ν^2*qr) (1+sqrt[4+(Zni/(4*qr))^2])])^0.5];
D0[ne_?NumericQ, qr_?NumericQ] := NIntegrate[16*Pi*νei*ne*Exp[-2*qr*r^2]*qr*r, {r, dmin, Infinity}];
A[ν_?NumericQ, qr_?NumericQ, ϵr_?NumericQ, qi_?NumericQ,ϵi_?NumericQ, dmin_?NumericQ] := 
  ν*NIntegrate[r*Exp[-qr*r^2]*((1-ϵr)*Sin[qi*r^2]-ϵi*Cos[qi*r^2]), {r, dmin, Infinity}];
B[ν_?NumericQ, qr_?NumericQ, ϵr_?NumericQ, qi_?NumericQ,ϵi_?NumericQ, dmin_?NumericQ] := 
  ν*NIntegrate[r*Exp[-qr*r^2]*((1-ϵr)*Cos[qi*r^2]+ϵi*Sin[qi*r^2]), {r, dmin, Infinity}];

s = NDSolve[{
  ν'[z] == -A[ν[z], qr[z], ϵr,qi[z], ϵi, dmin]*qr[z] -  D0[ne, qr[z]]*ν[z],
  qr'[z] == -2 A[ν[z], qr[z], ϵr,qi[z], ϵi, dmin]*qr[z]^2/ν[z]-D0[ne,qr[z]]*qr[z], 
  qi'[z] == -3 A[ν[z], qr[z], ϵr,qi[z], ϵi, dmin]*qr[z]*qi[z]/ν[z] +                        
               B[ν[z], qr[z], ϵr,qi[z], ϵi, dmin]*qr[z]^2/ν[z] -            
               D0[ne, qr[z]]*qi[z],
  ν[0.1] == 1., qr[0.1] == 10^(-4), qi[0.1] == 0.}, {ν, qr, qi}, {z, 0.1, 10000.}]
  • $\begingroup$ Your understanding for function is wrong. Just evaluate z[t, t, t] separately and think about why this happens. You may also want to read this question. Also, x[t] should be taken out of NIntegrate i.e. it should be x[t] NIntegrate[r*r, {r, 0, 5}]. $\endgroup$
    – xzczd
    Jun 13 '16 at 6:51
  • $\begingroup$ I am still trying to understand, but... And in my original code, x is a function of r, so it was x[r] not x[t]. $\endgroup$
    – sixpenny
    Jun 13 '16 at 7:11
  • $\begingroup$ Then I doubt if NDSolve will handle your code correctly. As far as I can tell, NDSolve isn't able to directly solve this type of integro-differential equation (at least currently). $\endgroup$
    – xzczd
    Jun 13 '16 at 7:27
  • $\begingroup$ If it is true, that is a big pity! I thought that mma is more powerful than matlab, then spent weeks to try to solve the integro-differential equations with mma. I still can not solve it. On the contrary, I only spent one day to solve the problem with matlab using runge-kuta method. And the coding with matlab is more straight forward. Maybe both mma and matlab have their own advantages. It will be useful to point out the weakness of mma, especially compared to matlab. Thus more guys will save their time! $\endgroup$
    – sixpenny
    Jun 13 '16 at 15:03
  • $\begingroup$ I began to use Mathematica because of a set of complicated PDEs, too, and as far as I can tell, for solving PDEs, Mathematica is much more handy then MATLAB. Though NDSolve can't directly handle this type of DE, Mathematica can solve several types of integro-differential equation (yeah there's a variety of integro-differential equations!) with just a little additional programming, for example this, this, this... $\endgroup$
    – xzczd
    Jun 13 '16 at 15:07

As mentioned in the comments above, your understanding for function is wrong. You'd better put more effort in learning the core language of Mathematica, or even consider revisiting the definition of function in math. The following is a not-the-simplest but correct way to code your equation set in Mathematica, do read it carefully:

vei = 1/1000;
r0 = Sqrt[1 + v^2 Exp[-2 qr rho^2]];
ne = 1/1000 - 4 qr r0 (1 - 1/r0^2) (1 - qr rho^2 (1 + 1/r0^2));
er = 1 - ne/r0; ei = 4 Pi vei ne;

d[v_, qr_] = Sqrt[(-Log@(1/1000 1/(v^2 qr) (1 + Sqrt[4 + (1/(1000 4 qr))^2])))]/Sqrt[qr];

{coreA[v_, qr_, qi_, rho_], coreB[v_, qr_, qi_, rho_]} = 
  Exp[-qr rho^2] rho RotationMatrix[qi rho^2].{ei, 1 - er};
int[core_, v_, qr_, qi_?NumericQ] := 
 v NIntegrate[core[v, qr, qi, rho] , {rho, d[v, qr], Infinity}, MaxRecursion -> 40]

coreD[v_, qr_, rho_] = 16 Pi vei ne Exp[-2 qr rho^2] qr rho;
intD[v_, qr_?NumericQ] := 
 NIntegrate[coreD[v, qr, rho], {rho, 0, Infinity}, MaxRecursion -> 40]

(* Notice I've actually modified the sign of A! *)    
eqn = With[{v = v@z, qr = qr@z, qi = qi@z}, 
   With[{intA = int[coreA, v, qr, qi], intB = int[coreB, v, qr, qi], 
     intD = intD[v, qr]}, {D[v, z] == intA qr - intD v, 
     D[qr, z] == 2 intA qr^2/v - intD qr, 
     D[qi, z] == 3 intA qi qr/v + intB qr^2/v - intD qi}]];

bc = {v[0] == 1, qr[0] == 10^(-4), qi[0] == 0};

{solv, solqr, solqi} = NDSolveValue[{eqn, bc}, {v, qr, qi}, {z, 0, 1}]
(* Some warning generated, but a result is still given *)
Plot[{Re@#[z], Im@#@z} &@solv // Evaluate, {z, 0, 1}]

enter image description here

When solving the equation set, NIntegrate spits out a lot of warning, I'm not sure about the reason. Perhaps it's a clue of divergence and the model should be checked, perhaps a higher WorkingPrecision or MaxRecursion is needed. Anyway I'd like to stop here and leave the remaining work to you.

  • $\begingroup$ I understood! I revised my code accordingly, and solved the problem. Inside the integration, core[v, qr, qi, rho] not 'core' should be used. While outside the integration, ne, not only ne[v, qr, qi, rho] can be used. This is one of the important issues I learned. Thank you very much! $\endgroup$
    – sixpenny
    Jun 22 '16 at 3:44
  • $\begingroup$ @sixpenny Actually the fatal part isn't the integration. What's crucial is, the independent variable must explictly exist in the body of function. Just compare the output of Clear[a, b, f]; b=a; f[a_]=b^2; f[2] and Clear[a, b, f]; b=a; f[a_]:=b^2; f[2] and think about why. If still feel confused, check the document ofSet and SetDelayed and try Clear[a, b, f]; b=a; Trace[f[a_]=b^2] $\endgroup$
    – xzczd
    Jun 22 '16 at 3:56

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.