ChrisRackauckas
commented
2 years ago noiseeqs = [Equation(p_noise, 1.0)]
That doesn't make sense. What are your actual noise equations?
Open RohanRajagopal opened 2 years ago
ChrisRackauckas
commented
2 years ago noiseeqs = [Equation(p_noise, 1.0)]
That doesn't make sense. What are your actual noise equations?
RohanRajagopal
commented
2 years ago Oh right, apologies. This is the stochastic system
$$\dot{x}_0(t) =x3(t),$$
$$\dot{x_1}(t) =x4(t),$$
$$\dot{x_2}(t)=x5(t),$$
$$\dot{x_3}(t) =Aa Sigm(x1(t)−x2(t))−2ax3(t)−a2x0(t),$$
$$\dot{x4}(t) =Aa[p{noise}(t)+C2Sigm(C1x0(t))]−2ax4(t)−a2x1(t),$$
$$\dot{x_5}(t) =BbC4Sigm(C3x0(t))−2bx5(t)−b2x2(t),$$
$$p_{noise} = p_0 + p_n + p_s sin(2 \pi p_s)$$
Here, $p_0$ is bias, $p_n$ is Gaussian white noise with zero mean and standard deviation, $\sigma_n$.
This is what I was trying to implement, but I think the $proper$ way to phrase it would have been the following(?),
eqs = [D(y0) ~ y3,
D(y1) ~ y4,
D(y2) ~ y5,
D(y3) ~ A * a * logistic(y_output) - 2 * a * y3 - a * a * y0,
D(y4) ~ A * a * (C2 * logistic(C1 * y0)) - 2 * a * y4 - a * a * y5,
D(y5) ~ B * b * C4 * logistic(C3 * y0) - 2 * b * y5 - b * b * y2,
y_output ~ y1 - y2
]
noiseeqs = [0,
0,
0,
0,
1.0,
0]Thank you.
(Added)
I believe, I have been able to implement it, at least in part correctly.
using ModelingToolkit, Symbolics, Plots, DifferentialEquations, LogExpFunctions
e_0 = 2.5
ν_0 = 6
r = 0.56
@parameters t, a, A, b, B, C_0, C1, C2, C3, C4, p_0, f_s, p_s
@variables y0(t), y1(t), y2(t), y3(t), y4(t), y5(t), y_output(t), p_noise3(t)
D = Differential(t)
Sigm(x)= @. 2*e_0 / (1 + exp(r*(ν_0 - x)))
@register Sigm(x)
p_noise1(t) = @. sinpi(p_0)^2 + sinpi(p_s)^2 * sinpi(2 * f_s * t)
p_noise2(t) = @. p_0 + p_s * sinpi(2 * f_s * t)
eqs = [D(y0) ~ y3,
D(y1) ~ y4,
D(y2) ~ y5,
D(y3) ~ A * a * Sigm(y1 - y2) - 2 * a * y3 - a * a * y0,
D(y4) ~ A * a * (p_noise2(t) + C2 * Sigm(C1 * y0)) - 2 * a * y4 - a * a * y1,
D(y5) ~ B * b * C4 * Sigm(C3 * y0) - 2 * b * y5 - b * b * y2,
]
noiseeqs = [0,
0,
0,
0,
32.5*4.8989,
0]
@named s_prob = SDESystem(eqs, noiseeqs, t, [y0, y1, y2, y3, y4, y5], [a, A, b, B, C1, C2, C3, C4, p_0, f_s, p_s])
#s_prob = structural_simplify(s_prob)
#@named D_prob = ODESystem(eqs, t, [y0, y1, y2, y3, y4, y5, y_output, p_noise3], [a, A, b, B, C1, C2, C3, C4, p_0, f_s, p_s])
#D_prob = structural_simplify(D_prob)
const c_0 = 132
p0 = [
a => 100
A => 3.25
b => 50
B => 22
C1 => c_0
C2 => c_0*0.8
C3 => c_0/4
C4 => c_0/4
p_0 => 125
f_s => 8.5
p_s => 20
]
y_0 = [
y0 => 5,
y1 => 5,
y2 => 5,
y3 => 5,
y4 => 5,
y5 => 5
]
tspan = (0.0, 400.0)
#prob_ode = ODEProblem(D_prob, y_0, tspan, p0)
dt = 0.0000001
prob_sde = SDEProblem(s_prob, y_0, tspan, p0)
sol = solve(prob_sde, LambaEulerHeun(), dt = dt)
plot(sol.t, sol[y1] - sol[y2])
When I use the adaptive solvers the $maxiters$ error is encountered with $dt$ getting smaller and smaller; is there anyway that can be taken care of. And, so I used the $LambaEulerHeun()$.
Thanks very much.
ChrisRackauckas
commented
2 years ago
I was asked to open an issue. Thank you very much.
Setting up a SDESystem, encountering an error while using structural_simplify that I can’t quite understand.
The error encountered is as follows
Thank you very much for your incredible work.