1
$\begingroup$

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.

Thanks in advance

$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – Lukas
    Commented Nov 15, 2015 at 19:54
  • $\begingroup$ 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,...}. $\endgroup$
    – holistic
    Commented Nov 15, 2015 at 20:17

1 Answer 1

1
$\begingroup$

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.

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

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.