I have to solve a series of differential equations along with a constraint that has to be satisfied at each point of time. To solve the differential equations, I first solves the following in order to obtain λ(t) and x(t)
sol1 = NDSolve[{λ'[t] == -((c1^2 - 4*k1*β*λ[t] -
2*c1*(λ[t] + μ) + (λ[t] + μ)^2)/(
4 k1)), x'[t] == -((c2^2 - 4*k2*(β*x[t] + δ1 (x[t] - [Lambda][t])) - 2*c2*(x[t] + μ) + (x[t] + μ)^2)/(4 k2)), λ[T] == 0, x[T] == 0}, {λ, x}, {t, 0, 100}];
Then using those values I calculate ρ1(t) and ρ2(t)
ρ1[t_] := (c1 - {λ[t] /. sol1} - μ)/(2*k1);
ρ1[t][[1, 1]];
ρ2[t_] := (c2 - {x[t] /. sol1} - μ)/(2*k2);
ρ2[t][[1, 1]];
Having ρ1(t) and λ(t), ρ2(t) and x(t) I solve the following differential equations where
d[t_] := a + b*Sin[(2*π*t)/25];
s = NDSolve[{y3'[t] == d[t] - δ2*y3[t], y2'[t] == δ2*y3[t] - (δ1 - ρ2[t][[1, 1]]) y2[t], y1'[t] == δ1*y2[t] - ρ1[t][[1, 1]]*y1[t],
y3[0] == Subscript[y, 0], y2[0] == Subscript[y, 0],
y1[0] == Subscript[y, 0]}, {y3, y2, y1}, {t, 0, 100}];
Now I have to check if the constraint is the satisfied, the constraint is the following: if
d[t] - (ρ1[t]*{y1[t] /. s} + ρ2[t]*{y2[t] /. s}) >= 0
then μ == 0
Otherwise calculate the μ(t) from below
μ == (-2*d[t]*k1*k2 + c1*k2*{y1[t] /. s} + c2*k1*{y2[t] /. s} - k2*{y1[t] /. s}*{λ[t] /. sol1} - k1*{y2[t] /. s}*{x[t] /. sol1})/(k2*{y1[t] /. s} + k1*{y2[t] /. s})]
and then I have to solve everything with the new value of μ.
The model parameters are
β = 0.1; c1 = 10; c2 = 30; k1 = 100; k2 = 300; T = 100; \
δ1 = 0.1; δ2 = 0.3; Subscript[y, 0] = 1000; a = 4000; b \
= 200;
Is there anyways to code this model to check for μ(t) and assign the value for it based on the defined condition? I also need to plot d[t] - (ρ1[t]*{y1[t] /. s} + ρ2[t]*{y2[t] /. s})
to see for what intervals my constraint is violated.
p.s. I'm using mathematica v.9.0.
Thanks in advanced for any help