issue with trapezoidal rule to solve IVP ODE, any lead appreciated

Problems with modeling
Post Reply
shoebmoon
User
User
Posts: 4
Joined: 5 years ago

issue with trapezoidal rule to solve IVP ODE, any lead appreciated

Post by shoebmoon »

Hi everyone,

Recently I am facing problem in modeling trapezoidal rule to solve initial value ODE through descretization .
Trapezoidal rule for IVP or BVP ODE :

dy/dx = f(y);
y(n+1) = y(n) + h/2 ( f(y(n+1)) + f(y(n))) /here h, is the step size/
y(1)=a


Example ODE: A*dh/dt = theta - kt*qout , where qout =sqrt(h) (so overall its a single variable ODE )
IV: h(1) = 1
Time Period: 800 min
Step size: 1 min


Below is the code that i wrote: sets

*Time horizon of 800 minutes
m /1*800/;

Scalar
A cross sectional area /5/
hinit initial height /1/
theta inlet flow rate /0.5/
kt constant /0.3/ ;

Variables
h(m)
hmax
qout(m) ;

Equations

htt
eq1(m)
eq2(m)
obj;

*Trapezoidal rule

eq1(m)$(ord(m) lt 800).. A*h(m+1)-h(m) =e= 1/2*(theta-qout(m+1)+theta-qout(m));

eq2(m)$(ord(m) lt 800).. qout(m)=e=kt*sqrt(h(m));

*Initial value
htt.. h('1')=e=hinit;

*Dummy objective
obj.. hmax =e= smax(m,h(m));

Model profile /all/ ;

*Variable bounds
h.lo(m)=1;
h.up(m)=5;
option dnlp=coinipopt;
solve profile minimizing hmax using dnlp;
display h.l; Also attached the code file below!

My requirement is to see the accuracy of approximation via Trapezoidal and for that purpose i want the profile of 'h' with time .
With the above code the value of dynamic variable 'h' is not changing from the initial value of 1, throughout the time interval.
I am making some semantic error.
Any help is much appreciated.

Thanks a lot!
Attachments
trapezoidiscrete.gms
(1.44 KiB) Downloaded 151 times
Post Reply