1
$\begingroup$

One example for coupled ODEs from the Mathemtica help center reads

a = NDSolve[{x'[t] == -3 (x[t] - y[t]), 
    y'[t] == -x[t] z[t] + 26.5 x[t] - y[t], 
    z'[t] == x[t] y[t] - z[t], x[0] == z[0] == 0, 
    y[0] == 1}, {x, y, z}, {t, 0, 5}]

Now we can for example extend this expression by replacing

x'[t] == -3 (x[t] - y[t]) -> x['t] == -3* (x[t] - y[t])*int[t]

with

f[t_] := Sin[t^2]
int[t_?NumericQ] := NIntegrate[Exp[-f[t]]*v^2, {v, 0, 1}]

This works indeed fine, but now I would like to extend the integral expression further by making the integrand dependent on the current solution of z[t] i.e.

int[t_?NumericQ] := NIntegrate[Exp[-f[t]*Evaluate[z[t]/.a]*v^2], {v, 0, 1}]

But this expression returns a bunch of errors. Does anyone know how to make this work?

$\endgroup$
9
  • $\begingroup$ But neither z nor f are dependent on v, no? $\endgroup$ Jul 25, 2016 at 13:41
  • $\begingroup$ Right, they are both independent $\endgroup$
    – stef
    Jul 25, 2016 at 13:59
  • $\begingroup$ Then you can factor them out of the integral, no? $\endgroup$ Jul 25, 2016 at 14:21
  • $\begingroup$ Thats not trivial if one has e.g. Sinh[f[t]*z[t]*v] as integrand $\endgroup$
    – stef
    Jul 25, 2016 at 14:55
  • 2
    $\begingroup$ You probably should have used that example, then. :) $\endgroup$ Jul 25, 2016 at 14:58

1 Answer 1

1
$\begingroup$

I think substituting a in the definition of int, which is called by a is trouble, and what you really need is to make sure int isn't evaluated when z[t] isn't a number. The following seems to work, albeit not very quickly:

Clear[a];

f[t_] := Sin[t^2];
int[t_?NumericQ, z_?NumericQ] := NIntegrate[Exp[-f[t]*z*v^2], {v, 0, 1}];

a = NDSolve[{
  x'[t] == -3*(x[t] - y[t])*int[t, z[t]], 
  y'[t] == -x[t] z[t] + 26.5 x[t] - y[t],
  z'[t] == x[t] y[t] - z[t], 
  x[0] == z[0] == 0, y[0] == 1},
  {x, y, z}, {t, 0, 5}]
$\endgroup$
2
  • $\begingroup$ That's what I was looking for, works fine! However, adding to the above definition of int[t,z] the new def. int2[t_?NumericQ, z_?NumericQ] := h[1] /. NDSolve[{h'[v] == Exp[-f[t]*z*v^2], h[0] == 0}, h, {v, 0, 1}] and calling int2[t,z[t]] in NDSolve produces fatal errors, even though Plot[{int[t, 1], int2[t, 1]}, {t, 0, 5}] shows the same result. Do you know how to do the replacement right? $\endgroup$
    – stef
    Jul 27, 2016 at 9:14
  • 1
    $\begingroup$ int2[2, 1] gives {1.32159} not 1.32159 -- you have to strip off a level with [[1]]. That is, int2[t_?NumericQ, z_?NumericQ] := h[1] /. NDSolve[{h'[v] == Exp[-f[t]*z*v^2], h[0] == 0}, h, {v, 0, 1}][[1]] works $\endgroup$
    – Chris K
    Jul 27, 2016 at 18:38

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.