Take the 2-minute tour ×

My questioin is: why rhsexpr works but rhsfunc doesn't. Thanks. :)

    m = l = 1.; \[Mu] = 100;
    beltv[t_] = .1;
    spring[x_] = 1000. (l - x);

    viscous[v_] := -30. (v - beltv[t]);
    coulomb[v_] := -25. Sign[v - beltv[t]];
    stribeck[v_] := -.3 Sign[v] Exp[-2 Abs[v]];

    friction[v_] := viscous[v] + coulomb[v] + stribeck[v];

rhsfunc is defined as:

    rhsfunc[t_, x_, 1] := beltv'[t];
    rhsfunc[t_, x_, 0] := spring[x[t]] + friction[x'[t]];

rhsexpr is defined as:

    rhsexpr := 
      If[stuck[t] == 1, beltv'[t], spring[x[t]] + friction[x'[t]]];

rhsexpr works.

     NDSolveValue[{x''[t] == rhsexpr , x[0] == 1, x'[0] == 0, 
     WhenEvent[x'[t] == beltv[t] && spring[x[t]]^2 <  \[Mu]^2, stuck[t] -> 1], 
      WhenEvent[spring[x[t]]^2 >=  \[Mu]^2, stuck[t] -> 0], 
      stuck[0] == 0}, {x[t], x'[t]}, {t, 0, 2}, 
     DiscreteVariables -> stuck[t]]

rhsfunc doesn't work, prompting that Encountered non-numerical value for a derivative at t == 0

    NDSolveValue[{x''[t] == rhsfunc[t, x[t], stuck[t]] , x[0] == 1, 
      x'[0] == 0, 
       stuck[t] -> 1], 
      WhenEvent[spring[x[t]]^2 >=  \[Mu]^2, stuck[t] -> 0], stuck[0] == 0}, {x[t], x'[t]}, {t, 0, 2},  DiscreteVariables -> stuck[t]]

enter image description here

update 1

Actually, I defines the rhsfunc just as the m@ does as in its bulit-in documentalion which is showed below:

rhs[t_, x_, 1] = -x;
rhs[t_, x_, 2] = 1;

sol = NDSolve[{x'[t] == rhs[t, x[t], a[t]], x[0] == 1, a[0] == 1, 
WhenEvent[x[t] - .1, a[t] -> 2], 
WhenEvent[x[t] - 2, a[t] -> 1]}, {x, a}, {t, 0, 15}, 
DiscreteVariables -> {a[t] \[Element] {1, 2}}];

Plot[x[t] /. sol, {t, 0, 10}]

enter image description here

share|improve this question
Try rhsfunc[t_?NumberQ, x_?NumberQ,1]:= ... –  Matariki Mar 24 '13 at 7:44
Thanks Matariki. I tried it, but it doesn't work, still prompting "non-numerical value for a derivative at t == 0". –  pengfei_guo Mar 24 '13 at 7:49
be sure to clear your old def... –  george Mar 24 '13 at 12:30
I am sure that I cleared everthing before I modify then run according to matariki's suggestion. I even quit the kernel before I run the modified code. Please see my attached screenshot. It just doesn't work. Anyway, thanks george. –  pengfei_guo Mar 24 '13 at 13:15
try rhsfunc[t_, x_, __] for the second case (instead of explicit 0). that is actually what your rhsexpr does. I'm guessing on first pass ndsolve tries a symbolic value for stuck[] which is undefined as you have it. –  george Mar 24 '13 at 19:29
show 5 more comments

1 Answer

up vote 1 down vote accepted

Patterns are not much like If statement, especially in your example, where the If ststement is delayed. Patterns require an exact agreement with the definitions, and since you didn't define what happens for rhsfunc when it gets non numerical value for the third argument, it remain unevaluated and NDSolve cannot proceed its evaluation The last example provided differs from yours since it does not involves derivatives, and in addition, the call is different in format (e.g., returning -x rather -x[t] and when you call it with x[t] it return the pattern -x[t])

One additional note - NDSolve evaluates its argument before it functions on them (since it has no HoldAll or HoldFirs or HoldRest attribute). This means that the rhsfunc is evaluated immediately with a symbol returning a non evaluated expression As a rule of thumb, when using NDSolve prefer always Piecewise/ If/ Which. Piecewise was developed to function together with NDSolve and alike functions because of this behavior (starting in ver. 5.1). If you use patterns, use ?NumericQ as a test for the variables (NumberQ does not return True for expressions involving Pi or E etc., while NumericQ does). This prevents the function from being evaluated for non numeric variables

I hope this helps


share|improve this answer
add comment

Your Answer


