Prev: How can I plot two curves with two different y-axis (placed at right and left) with one reverse-logarithmic x-axis?
Next: array of struct like dir function
From: Floris Zoutman on 29 Apr 2010 08:51
I am solving a system of dynamic and algebraic equations using ode23t in combination with a (singular) mass matrix. One of my variables has a non-negativity constraint and I have some problems programming this. In mathematical form my equations look like this:
dy/dt=f(t,x,y,p) for all t
0=g(t,x,y,p,eta) for all t
0=p*eta for all t
p(t)>=0 for all t
eta(t)>=0 for all t
x(0)=x0
y(0)=y0
p(0)=p0>0
eta(0)=0
x(t) is a n*1 vector, y is a m*1 vector, p is the (scalar) variable that has a non-negativity constraint and eta its Kuhn-Tucker multiplier, f maps to R^m and g maps to R^(n+1) (one extra for p).

If you exclude the inequality constraints there are two time-paths that solve the equations. In the first p is strictly decreasing and eta =0 for all t. In this case p becomes negative at the end-point. This is not an interesting time-path for me because p cannot be negative. The other solution is one in which p is strictly decreasing in t until it hits zero. After that p stays at zero and eta becomes positive. This is the time path I am interested in. Unfortunately, I do not see how to force Matlab to solve for the second instead of the first time path. Therefore my question is: can you program ode23t to honour my non-negativity constraints? Alternatively, is there a smart mathematical trick to circumvent this problem?
From: Torsten Hennig on 29 Apr 2010 22:30
> I am solving a system of dynamic and algebraic
> equations using ode23t in combination with a
> (singular) mass matrix. One of my variables has a
> non-negativity constraint and I have some problems
> programming this. In mathematical form my equations
> look like this:
> dy/dt=f(t,x,y,p) for all t
> 0=g(t,x,y,p,eta) for all t
> 0=p*eta for all t
> p(t)>=0 for all t
> eta(t)>=0 for all t
> x(0)=x0
> y(0)=y0
> p(0)=p0>0
> eta(0)=0
> x(t) is a n*1 vector, y is a m*1 vector, p is the
> (scalar) variable that has a non-negativity
> constraint and eta its Kuhn-Tucker multiplier, f maps
> to R^m and g maps to R^(n+1) (one extra for p).
>
> If you exclude the inequality constraints there are
> two time-paths that solve the equations. In the
> first p is strictly decreasing and eta =0 for all t.
> In this case p becomes negative at the end-point.
> This is not an interesting time-path for me because p
> cannot be negative. The other solution is one in
> which p is strictly decreasing in t until it hits
> zero. After that p stays at zero and eta becomes
> positive. This is the time path I am interested in.
> Unfortunately, I do not see how to force Matlab to
> solve for the second instead of the first time path.
> Therefore my question is: can you program ode23t to
> honour my non-negativity constraints? Alternatively,
> is there a smart mathematical trick to circumvent
> this problem?

Substitute
eta(t) := eta~(t)^2, p(t):=p~(t)^2
(or use some other non-negative function for
substitution).

Best wishes
Torsten.
 | 
Pages: 1
Prev: How can I plot two curves with two different y-axis (placed at right and left) with one reverse-logarithmic x-axis?
Next: array of struct like dir function