# Differential equation for a list with parameter dependent function

I am having a differential equation:

y' = (1 - y) - f[y, mu] y;


f is a hysteretic function that depends on y and on the derivative of y:

f[y_, mu_] := -1/(1 + Exp[(Abs[y] - mu)/k]) + 1;
mu[n_] := mu[n] = 0.5 + 0.1 Sign[y[n] - y[n - 1]];


For y being a scalar, I can solve the equation numerically:

x[0] = 0.0; y[0] = 0; y[1] = 0.1; x[n] = 2; h = .01;
x[n_] := x[n] = x[0] + n h;

y[n_] := Module[{k1, k2, k3, k4},
k1 = h (f[x[n - 1], y[n - 1], mu[n - 1]]);
k2 = h (f[x[n - 1] + h/2, y[n - 1] + k1/2, mu[n - 1]]);
k3 = h (f[x[n - 1] + h/2, y[n - 1] + k2/2, mu[n - 1]]);
k4 = h (f[x[n - 1] + h, y[n - 1] + k3, mu[n - 1]]);
y[n] = y[n - 1] + (k1 + 2 k2 + 2 k3 + k4)/6];


However, I need to solve the equation for a vector y with y[0]={0,0}. How can I expand Runge - Kutta for a list and still include the hysteretic function?

• I love that you made your own RK4 algorithm, it's so much more fun than using the ExplicitRungeKutta method for NDSolve. But here, I don't see why you should have any problem using a list for y. If you take your f function and give it a list for the y argument, and a list for the mu argument, it returns a list, since both Exp and Abs thread over lists. Jan 19, 2016 at 11:04
• But I do see a problem with the code in that you are always passing 3 arguments to f, when you have f defined as taking only 2 arguments. Also, k is undefined Jan 19, 2016 at 11:07
• Greetings! To make the most of Mma.SE please take the tour now. Help us to help you, write an excellent question. Edit if improvable, show due diligence, give brief context, include minimum working examples of code and data in formatted form. As you receive give back, vote and answer questions, keep the site useful, be kind, correct mistakes and share what you have learned. Jan 19, 2016 at 11:39
• You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is also useful for learning how to format your questions and answers. You may also find this this meta Q&A helpful Jan 19, 2016 at 12:14

I don't know what value of k you are using, so in the code below I set it to 0.1. Also, fixed the RK4 code so that it works, but now it's pretty slow.

(* set k to your value for k *)
k = .1;
f[y_, mu_] := -1/(1 + Exp[(Abs[y] - mu)/k]) + 1;
mu[n_] := mu[n] = 0.5 + 0.1 Sign[y[n] - y[n - 1]]; x[0] = 0.0;
y[0] = 0; y[1] = 0.1; x[n] = 2; h = .01;
x[n_] := x[n] =
x[0] + n h;(*Why do you have x defined? Neither f nor mu depend on \
x*)

y[n_] := Module[{k1, k2, k3, k4},
k1 = h (f[y[n - 1], mu[n - 1]]);
k2 = h (f[y[n - 1] + k1/2, mu[n - 1]]);
k3 = h (f[y[n - 1] + k2/2, mu[n - 1]]);
k4 = h (f[y[n - 1] + k3, mu[n - 1]]);
y[n - 1] + (k1 + 2 k2 + 2 k3 + k4)/6];

Do[
Print[AbsoluteTiming@y[n]]
, {n, 1, 10}]

(*
{2.*10^-6,0.1}
{0.000128,0.100067}
{0.000468,0.100134}
{0.00234,0.100201}
{0.010894,0.100268}
{0.038722,0.100335}
{0.145408,0.100402}
{0.70548,0.10047}
{3.46992,0.100537}
{17.5871,0.100604}
*)


The time to compute the next y value is increasing exponentially. The answer here is to make y a function of yin:

k = .1;
f[y_, mu_] := -1/(1 + Exp[(Abs[y] - mu)/k]) + 1;
ycomp = Compile[{{yin, _Real}, {muin, _Real}, {step, _Real}},
Module[{k1, k2, k3, k4, yout, muout},
k1 = step (f[yin, muin]);
k2 = step (f[yin + k1/2, muin]);
k3 = step (f[yin + k2/2, muin]);
k4 = step (f[yin + k3, muin]);
yout = yin + (k1 + 2 k2 + 2 k3 + k4)/6;
muout = 0.5 + 0.1 Sign[yout - yin];
{yout, muout}]]


Now you get the exact same values for y in less than a millisecond.

{yc, muc} = {0.1, 0.6};
AbsoluteTiming@Table[{yc, muc} = ycomp[yc, muc, .01]; yc, {10}]
(* {0.000194, {0.100067, 0.100134, 0.100201, 0.100268, 0.100335,
0.100402, 0.10047, 0.100537, 0.100604, 0.100672}} *)


So now, what if y were a list? Just modify the Compile code a little bit,

k = .1;
f[y_, mu_] := -1/(1 + Exp[(Abs[y] - mu)/k]) + 1;
ycomp = Compile[{{yin, _Real, 1}, {muin, _Real, 1}, {step, _Real}},
Module[{k1, k2, k3, k4, yout, muout},
k1 = step (f[yin, muin]);
k2 = step (f[yin + k1/2, muin]);
k3 = step (f[yin + k2/2, muin]);
k4 = step (f[yin + k3, muin]);
yout = yin + (k1 + 2 k2 + 2 k3 + k4)/6;
muout = 0.5 + 0.1 Sign[yout - yin];
{yout, muout}]]


Now all you have to do is supply the initial yc and muc as lists of any length and all is good.

yc = {0.1, 0.11};
muc = {0.6, 0.66};
AbsoluteTiming@Table[{yc, muc} = ycomp[yc, muc, .01]; yc, {10}]
(* {0.000895, {{0.100067, 0.110041}, {0.100134,
0.110115}, {0.100201, 0.110189}, {0.100268, 0.110263}, {0.100335,
0.110337}, {0.100402, 0.110411}, {0.10047, 0.110485}, {0.100537,
0.11056}, {0.100604, 0.110634}, {0.100672, 0.110708}}} *)

• thank you very much. The code works very well. Jan 20, 2016 at 14:09
• Great, glad I could help. If that solves the problem, go ahead and mark it as solved so it doesn't stay in the "Unanswered questions" queue. Jan 20, 2016 at 14:10