2
$\begingroup$

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]rr, {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_] := xr; v[u_?NumericQ] := NIntegrate[ur, {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_] := xr; v[x_?NumberQ] := NIntegrate[ur, {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_] := xr; v[x_?NumberQ] := NIntegrate[ur, {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 variables 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.

Equations:

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.}]
Improve this question
$\endgroup$
20
  • $\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$ Commented Jun 13, 2016 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$ Commented Jun 13, 2016 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$ Commented Jun 13, 2016 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$ Commented Jun 13, 2016 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$ Commented Jun 13, 2016 at 15:07

1 Answer 1

Reset to default
2
$\begingroup$

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:

Clear["`*"] 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.

Improve this answer
$\endgroup$
2
  • $\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$ Commented Jun 22, 2016 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$ Commented Jun 22, 2016 at 3:56

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.