7
$\begingroup$

I am trying to solve a linear ODE with a variable coefficient which is given in terms of an integral I can only do numerically. That is, I have an equation of the form $$ u'(t)=u(t)\int_0^t\sqrt{1+e^{-\tau^2}\sin(\tau)}d\tau, $$ where the integrand could be pretty much any well-behaved function - but which in general will cause Mathematica to choke when given to Integrate and eventually say it can't do it.

I know that such a lack of an analytic handle on the coefficients in the equation essentially disqualifies me from an analytic solution of the equation. However, given that I am, or should be, perfectly able to numerically compute this coefficient, Mathematica ought to be able to numerically solve the equation.

My initial guess as to how to do this would be to simply indicate the integral should be done numerically:

NDSolve[{D[u[t], t] == NIntegrate[
     Sqrt[1 + E^tt^2 Sin[tt^2]], {tt, 0, t}] u[t], u[0] == 1}, u, {t, 0, 1}]

but this simply outputs the error

NIntegrate::nlim: tt = t is not a valid limit of integration. >>.

I've also tried with plain Integrate and setting GenerateConditions->False, but none of them work and they also take a long time to return their errors, which I assume is because they're trying to symbolically solve the integral - and I specifically don't want that.

Can someone provide a way to make Mathematica understand it just needs to buckle down and calculate?

$\endgroup$
7
  • 2
    $\begingroup$ Define the integral as a separate function that can only evaluate for numerical arguments. For instance, f[t_?NumericQ] := NIntegrate[Sqrt[1 + Exp[-tt^2] Sin[tt^2]], {tt, 0, t}]; NDSolve[{u'[t] == f[t] u[t], u[0] == 1}, u, {t, 0, 1}] $\endgroup$ Commented Oct 26, 2012 at 16:56
  • $\begingroup$ A bracket was missing, edited. $\endgroup$ Commented Oct 26, 2012 at 17:02
  • $\begingroup$ @J.M. this question/issue really pops up every 2 weeks! $\endgroup$
    – chris
    Commented Oct 26, 2012 at 17:09
  • $\begingroup$ related (amongst many others) mathematica.stackexchange.com/a/13048/1089 $\endgroup$
    – chris
    Commented Oct 26, 2012 at 17:13
  • $\begingroup$ @episanty Why do you want to use tt to describe \[Tau]. Hit escape key , t and then escape again to get tau. I see this tt a lot with matlab and C users alike but this isn't needed in mathematica. $\endgroup$
    – dearN
    Commented Oct 26, 2012 at 17:26

3 Answers 3

7
$\begingroup$

Try

ClearAll[f];
f[t_Real] := NIntegrate[Sqrt[1 + E^tt^2 Sin[tt^2]], {tt, 0, t}];
NDSolve[{D[u[t], t] \[Equal] f[t]* u[t], u[0] == 1}, u, {t, 0, 1}]

ie, define f such that it only evaluates its right hand side if its argument has head Real.

$\endgroup$
10
$\begingroup$

As an alternative to the answer of @acl, let's approach the ODE from mathematical side.

Approach 1
Let's denote $v(t) = \int_0^t f(\tau) \mathrm{d} \tau$,where $f(\tau) = \sqrt{1 + \mathrm{e}^{-\tau^2} \sin(\tau)}$. Then clearly $v^\prime(t) = f(t)$ and $v(0) = 0$. You can thus rewrite you ODE as a system: $$ u^\prime(t) = v(t) u(t), \, v^\prime(t) = f(t), \, u(0) = 1, \, v(0) = 0 \tag{1} $$

Approach 2
Rewrite your equation as $$ \frac{\mathrm{d}}{\mathrm{d} t} \log u(t) = \frac{u^\prime(t)}{u(t)} = \int_0^t f(\tau) \mathrm{d}\tau \tag{2} $$ thus $$ u(t) = u(0) \exp\left( \int_0^t \left( \int_0^{s} f(\tau) \mathrm{d} \tau\right) \mathrm{d} s\right) = u(0) \exp(w(t)) \tag{exact} $$ where $w(t)$ satisfies $w^{\prime\prime}(t) = f(t)$ with initial conditions $w(0) = w^\prime(0) = 0$.

Approach 3
Alternatively, we could differentiate eq. (2) once more, and get $$ \frac{u^{\prime\prime}(t) u(t) - u^\prime(t) u^\prime(t)}{u(t)^2} = f(t), \quad \text{or} \quad u^{\prime\prime}(t) u(t) = f(t) u(t)^2 + \left( u^\prime(t)\right)^2 \tag{3} $$ with initial conditions $u^{\prime}(0) = 0$.

Let's compare three approaches using Mathematica:

enter image description here


Here are the inputs:

f[t_] := Sqrt[1 + Exp[-t^2] Sin[t]]

sol1 = NDSolve[{u'[t] == v[t] u[t], v'[t] == f[t], u[0] == 1, 
   v[0] == 0}, u, {t, 0, 4 Pi}]

sol2 = NDSolve[{u[t] u''[t] == f[t] u[t]^2 + u'[t]^2, u[0] == 1, 
   u'[0] == 0}, u, {t, 0, 4 Pi}]

solaux = NDSolve[{w''[t] == f[t], w[0] == 0, w'[0] == 0}, 
  w, {t, 0, 4 Pi}]

uexact[t_Real] := Exp[w[t] /. solaux]

LogPlot[{u[t] /. sol1, u[t] /. sol2, uexact[t]} // Evaluate, {t, 0, 
  4 Pi}, PlotStyle -> {Directive[Orange, Thickness[0.02]], 
   Directive[Black, Thickness[0.01]], Directive[Yellow, Dashed]}]
$\endgroup$
3
  • $\begingroup$ I appreciate the approach, but in my particular case if it's a re-think on the maths side then it needs to go deeper (which, for efficiency reasons, you need to do anyway). My specific problem was getting it to evaluate when needed and not complain when not, which is solved by the t_Real trick. $\endgroup$ Commented Nov 5, 2012 at 14:19
  • $\begingroup$ @episanty Notice that approach 1 does exactly that, but avoid calling NIntegrate, thus achieving much better computational complexity. $\endgroup$
    – Sasha
    Commented Nov 5, 2012 at 15:23
  • $\begingroup$ @Sasha, do you have a timing for those three methods? I'm looking at a problem which contains three nasty integrals in my differential equation. Is there any way to know which method will be the fastest ? $\endgroup$
    – Nick
    Commented Aug 13, 2015 at 15:41
2
$\begingroup$

First solve this:

a = NIntegrate[Sqrt[1 + E^tt^2 Sin[tt^2]], {tt, 0, 1}]

And then do this:

DSolve[{D[u[t], t] == a u[t], u[0] == 1}, u[t], t]

That should solve the problem.

$\endgroup$
4
  • $\begingroup$ the integral of the OP is indefinite: i.e. a=a[t]=Integrate[Sqrt[1 + E^tt^2 Sin[tt^2]], {tt, 0, t}] $\endgroup$
    – chris
    Commented Oct 26, 2012 at 18:34
  • $\begingroup$ @chris Oh crap! After doing this again, I came up with the same thing as @JM. But thanks for the -1 ;-) $\endgroup$
    – dearN
    Commented Oct 26, 2012 at 19:58
  • $\begingroup$ sorry I wanted to remove the -1 but its not possible; I just wanted to warn readers this answer fell short of answering the OP's question. $\endgroup$
    – chris
    Commented Oct 26, 2012 at 20:10
  • $\begingroup$ @chris Thats alright! Theres more to life than forum points! :P $\endgroup$
    – dearN
    Commented Oct 26, 2012 at 21:51

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.