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!