Question about Simpson's integration.

GAMS
1

Hello everyone. I have recently made a model to maximize the distance traveled by a glider. In this program, I have used the trapezoidal rule. However, my proffesor thinks that I can do it with Simpson’s rule instead. I have been given a small .pdf which can be sumarized (using LATEX code) as:

y_{k+1}= Y_{k}+\frac{h_{k}}{6}[f_{k+1}+4\hat{f_{k+1}}+f_{k}], where \hat{f_{k+1}=f(\hat{y_{k+1}}, \hat{u_{k+1}}, t_{k}+\frac{h_k}{2})} and \frac{y_{k+1}}= \frac{h_{k}}{2}(y_{k+1}+y_{k})+\frac{k_{k}{8}}(f_{k}-f{k+1})

And using regular script:

yk+1 = yk + hk/6[fk+1+4fk+1+fk], where fk+1=f(yk+1, uk+1, tk+hk/2) and yk+1= 1/2(yk+1+ yk)+hk/8(fk-fk+1).

The program I have written is stated next:

$set n 50
set j /0*%n%/;
sets
jlast(j)
jnotlast(j);
jlast(j)$(ord(j)=card(j))=yes;
jnotlast(j)=not jlast(j);
scalar
n number of intervals /%n%/
m mass /5000/
S surface /21.55/
CD0 drag /0.023/
k no idea /0.073/
hmax initial height /1000/
g gravity /9.81/

density density /1.225/
variable
gamma(j),
CL(j),
D(j),
CD(j),
L(j),
*x(j),
*y(j),
objective;

positive variable

x(j),
y(j),
v(j),
equation
diffx(j),
diffy(j),
valueD(j),
valueL(j),
obj;


 diffx[j]$(jnotlast(j)).. x[j+1]-x[j] =e=0.5*step*(v(j+1)*cos(gamma(j+1)) +  v(j)*cos(gamma(j))   );
 diffy[j]$(jnotlast(j)).. y[j+1]-y[j] =e=0.5*step*(v(j+1)*sin(gamma(j+1)) +  v(j)*sin(gamma(j))   );

valueD[j].. m*g*sin(gamma(j))=e=0.5*density*S*v(j)*v(j)*(CD0+k*CL(j)*CL(j));
valueL[j].. m*g*cos(gamma(j))=e=0.5*density*S*v(j)*v(j)*CL(j);

obj .. objective =e=   x('%n%');

x.fx('0') = 1.0e-12;
y.fx('0') = 1000;
y.fx('%n%') = 1.0e-12;

CL.up(j) =1.4;

y.up (j) = 1000;

gamma.up(j) = pi*0.5;

v.lo(j) =  1.0e-12;

y.lo(j) = 1.0e-12;

CL.lo(j) = 0;
gamma.lo(j) = 0;

model brahstron1 /all/;
nlp=ipopt;
solve brahstron1 using nlp maximize objective;

My issue is that since I don’t know how long it takes to the glider to reach it’s endpoint, I leave step as a free variable and I don’t know how to introduce the midpoints necessary to implement Simpson’s rule.

Any advice on how to proceed is welcome.
Thanks for reading.

2

If you find that using Simpson’s rule is throwing you off, then by all means implement a solution using the trapezoid rule. That will give you something, and you can move from that to Simpson’s quite easily, I believe. It seems to me the real difficulty is independent of which integration method is used.

Without seeing the homework assignment as it is given in the PDF, it’s difficult to say more.

-Steve

3

Edit: I have already found what was wrong. It was a licence problem; the functioning code is as follows:

$set n 10
set j /0*%n%/;
sets
jlast(j)
jnotlast(j);
jlast(j)$(ord(j)=card(j))=yes;
jnotlast(j)=not jlast(j);
scalar
n number of intervals /%n%/
m mass /5000/
S surface /21.55/
CD0 drag /0.023/
k ni idea /0.073/
hmax initial height /1000/
g gravity /9.81/

density density /1.225/
variable
gamma(j),
CL(j),
D(j),
CD(j),
L(j),

gamma_med(j),
CL_med(j),
D_med(j),
CD_med(j),
L_med(j),

objective;

positive variable

x(j),
y(j),
v(j),

x_med(j),
y_med(j),
v_med(j),

step;

equation
diffx(j),
diffy(j),
diffx_central(j),
diffy_central(j),
valueD(j),
valueL(j),
valueD_central(j),
valueL_central(j),
obj;


diffx[j]$(jnotlast(j)).. x[j+1]-x[j] =e=(1/6)*step*(v(j+1)*cos(gamma(j+1)) +  v(j)*cos(gamma(j)) + 4*v_med(j+1)*cos(gamma_med(j+1))  );
diffy[j]$(jnotlast(j)).. y[j+1]-y[j] =e=(-1)* (1/6)*step*(v(j+1)*sin(gamma(j+1)) +  v(j)*sin(gamma(j)) +  4*v_med(j+1)*sin(gamma_med(j+1))  );

diffx_central[j]$(jnotlast(j)).. x_med[j+1] =e=0.5*(x(j+1)+x(j)+(step/8)*(v_med(j)*cos(gamma_med(j)))-(v_med(j+1)*cos(gamma_med(j+1))));
diffy_central[j]$(jnotlast(j)).. y_med[j+1] =e=0.5*(y(j+1)+y(j)+(step/8)*(v_med(j)*sin(gamma_med(j)))-(v_med(j+1)*sin(gamma_med(j+1))));


valueD[j].. m*g*sin(gamma(j))=e=0.5*density*S*v(j)*v(j)*(CD0+k*CL(j)*CL(j));
valueL[j].. m*g*cos(gamma(j))=e=0.5*density*S*v(j)*v(j)*CL(j);

valueD_central[j].. m*g*sin(gamma_med(j))=e=0.5*density*S*v_med(j)*v_med(j)*(CD0+k*CL_med(j)*CL_med(j));
valueL_central[j].. m*g*cos(gamma_med(j))=e=0.5*density*S*v_med(j)*v_med(j)*CL_med(j);

obj .. objective =e=   x('%n%');




x.fx('0') = 1.0e-12;
y.fx('0') = 1000;
y.fx('%n%') = 1.0e-12;

CL.up(j) =1.4;
CL_med.up(j) =1.4;

y.up (j) = 1000;
y_med.up (j) = 1000;

gamma.up(j) = pi*0.5;
gamma_med.up(j) = pi*0.5;

gamma.lo(j) = 0;
gamma_med.lo(j) = 0;

v.lo(j) =  1.0e-12;
v_med.lo(j) =  1.0e-12;

y.lo(j) = 1.0e-12;
y_med.lo(j) = 1.0e-12;


CL.lo(j) = 0;
CL_med.lo(j) =0;

gamma.lo(j) = 0;
gamma_med.lo(j) = 0;


model brahstron1 /all/;
* Invoke the LGO solver option for solving this nonlinear programming
option
nlp=ipopt;
solve brahstron1 using nlp maximize objective;

I have just reduced the number of steps from 50 to 10.