# How to substitute a random function into a solution of a differential equation?

I would like to solve a differential equation including a source term.

The source term is a random variable and I model this by creating a list of RandomReal and then construct a 1-dim Interpolation to get a piecewise linear function.

I then try to replace all the constants appearing in the differential equation as well as the generic input function by values and my linear function.

Unfortunately, this results in a function that includes an Inactivated Integral of the Interpolation. I don't know, why I can't activate it and wether there is a deeper root of the error, that I don't understand.

randomValues = RandomReal[1, 100];
randomInput := Interpolation[randomValues, InterpolationOrder -> 1];
Plot[randomInput[x], {x, 1, 100},
PlotLabel -> "Piecewise linear function of random input"]
Print["Trying to solve the differential equation:"]
{U} = {u} /. DSolve[{u'[t] == u[t] + J[t], u[0] == 0}, {u}, t] //
FullSimplify // First
Print["Replacing constants and the piecewise Input Function"]
U = U /. {c1 -> 1., c2 -> 0., C -> 3., R -> 4., J -> randomInput}
Print["There is an inactive Integral involved. Trying to activate it."]
U = Activate[U]

U /@ Range[1, 100] (*List of Unevaluated Functions*)
Plot[U, {x, 1, 100}] (*No Error, but I am suspicious, wether it is the right output now*)
$$$$


I think it is better to use NDSolve for this. Having to integrate random function is not practical. (after activate, etc...) and will be very slow.

I used the noise[t] function from this answer Continuous noise representation and plugged it as input to the ode. Here is second order ode example using it (you can change the ode to the first order one you had ofcourse).

ClearAll[sol, u, t];
tmax = 10;
sigma = 1;
noise = Interpolation[
Normal[RandomFunction[WhiteNoiseProcess[sigma], {0, tmax}]][[1]]];
Plot[noise[t], {t, 0, tmax}, PlotLabel -> "Input, white noise"]


Now use NDSolve

sol = NDSolve[{u''[t] + u'[t] + u[t] == noise[t], u[0] == 2,
u'[0] == 0}, u, {t, 0, tmax}];

Plot[Evaluate[u[t] /. sol], {t, 0, tmax}, PlotLabel -> "Solution to ode"]


You can change the noise[t] to meet your noise specification. This will also work with DSolve but I found it very very slow to integrate.

To use DSolve

sol = Activate[u[t] /. First@ DSolve[{u''[t] + u'[t] + u[t] ==
noise[t], u[0] == 2, u'[0] == 0}, u[t], t]]


Plot[sol, {t, 0, tmax}]


It is very slow to plot compared to NDSolve plotting since of all the integrals there but if you wait, it will finish.

• Thank you. Nice work-around that also improved the evaluation time. But what was the core problem with the inactive integral? And how would I solve it? Nov 12 '21 at 12:59
• @Uwe.Schneider You just need to Activate` that is all. But it is slow. Added the code above. The slow part is the plotting, as it has to evaluate the intergals with that noise[t] function in there. Nov 12 '21 at 13:04