# Solving a differential equation with a specified "input vector"

I have the following differential equation:

sol = NDSolve[{v'[t] == (1)*(f[t] - v[t]^(1/0.2)),
q'[t] ==
f[t]*((1 - (1 - 0.8)^(1/f[t]))/0.8) - (v[t]^(1/0.2)*(q[t]/v[t])),
s'[t] == 0.5 u[t, 1, 2, 3, 4] - (s[t]/0.8) - (f[t] - 1)/0.4,
f'[t] == s[t], s[0] == 0, q[0] == v[0] == f[0] == 1}, {q, v, s,
f}, {t, 0, 30}]

where u[t,1,2,3,4] is a trapezoidal "stimulus" function serving as input to s[t]:

u[t_, a_, b_, c_, d_] :=
Piecewise[{{0, t < a}, {(t - a)/(b - a), a <= t <= b}, {1,
b <= t <= c}, {(d - t)/(d - c), c <= t <= d}, {0, d <= t}}]

So instead of u being a trapezoidal function, I want it to be some vector I specify. So I tried to replace u[t] with:

u={0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0,
0, 1, 0, 0, 0, 1, 1}

and used it in the following way, which doesn't work:

sol = NDSolve[{v'[t] == (1)*(f[t] - v[t]^(1/0.2)),
q'[t] ==
f[t]*((1 - (1 - 0.8)^(1/f[t]))/0.8) - (v[t]^(1/0.2)*(q[t]/v[t])),
s'[t] == 0.5 u[[t]] - (s[t]/0.8) - (f[t] - 1)/0.4,
f'[t] == s[t], s[0] == 0, q[0] == v[0] == f[0] == 1}, {q, v, s,
f}, {t, 0, 30}]

I get the following error message:

Part::pkspec1: The expression t cannot be used as a part specification. >>

Part::pkspec1: The expression t cannot be used as a part specification. >>

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`. >>

Can anybody give me some hints how to make this work? I'm sure it is possible since I saw an implementation of this equation in Matlab which used a prespecified vector rather than first defining a piecewise function.

• Since you do not specify what the problem actually is, my guess is the following: t will not be limited to integer values. So u[[t]] will cause trouble if evaluated for non-integers. You could probably define a function u that does always return a random number. Commented Nov 15, 2015 at 19:54
• Defining u[t_]:=RandomReal[{0,1}] works very well, thank you! I edited the thread because I not only want it to be a random number but also some vector I specify, e.g. u={u={0,0,1,0,1,0,1,1,...}. Commented Nov 15, 2015 at 20:17

It looks like your problem is that you are trying to define u as a list instead of as a function. You can turn your list into a function using interpolation. For your case:

uVals = {0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1,
0, 1, 0, 0, 1, 0, 0, 0, 1, 1};
u = Interpolation[Transpose[{Range[0, Length[uVals] - 1], uVals}]];
ϵ = 0.1;
sol = NDSolve[{v'[t] == (1)*(f[t] - v[t]^(1/0.2)),
q'[t] == f[t]*((1 - (1 - 0.8)^(1/f[t]))/0.8) - (v[t]^(1/0.2)*(q[t]/v[t])),
s'[t] == ϵ u[t] - (s[t]/0.8) - (f[t] - 1)/0.4,
f'[t] == s[t], s[0] == 0, q[0] == v[0] == f[0] == 1}, {q, v, s, f}, {t, 0, 29}]

It was also necessary to define epsilon, for which I arbitrarily chose the value 0.1.

• That worked great, thank you very much! I forgot to specify epsilon, but edited the post now. 0.1 was close enough though :). Commented Nov 16, 2015 at 8:53