Message Boards Message Boards

How to solve a constraint deferential equations with If function


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

POSTED BY: afshar_saman
2 years ago

You do not give the initial value used for μ in your first NDSolve. I assume there was some trial value and I picked zero.

You seem to be using { and } in an unconventional way to try to limit the scope of variable substitution. { and } are almost exclusively used for list construction. You might try using ( and ) instead.

For what seems to be the essential bit of your question, you might look at using

FindMinimum[d[t] - (ρ1[t]*(y1[t] /. s) + ρ2[t]*(y2[t] /. s)),{t, 50, 0, 100}]

which should give you a non-negative value if your constraint is satisfied and a negative value otherwise. Look at the information hidden behind Details and Options on the documentation page for FindMinimum to see how this should restrict the search for the minimum to lie between 0 and 100. That is only guaranteed to find local minima. There are a few other functions in Mathmatica, but I don't think they may guarantee a global minima with several interpolating functions. A function that has an option of using random search might stumble onto a minima less than zero when others might be trapped in local minima greater than zero.

I was able to plot the function you desired by using

Plot[First[d[t] - (\[Rho]1[t]*(y1[t] /. s) + \[Rho]2[t]*(y2[t] /. s))], {t, 0, 100}]

but I can't be certain that the use of First has not missed essential behavior perhaps hiding behind other solutions.

Hopefully you will see how you might adapt this to serve your purposes. All this should work with version 9.0.

POSTED BY: Bill Simpson
2 years ago

Group Abstract Group Abstract