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:
Code: Select all
*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;
Code: Select all
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!