Question about Simpson's integration.

Problems with modeling
Post Reply
confused student
User
User
Posts: 3
Joined: 5 years ago

Question about Simpson's integration.

Post by confused student »

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:

Code: Select all

$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.
User avatar
dirkse
Moderator
Moderator
Posts: 214
Joined: 7 years ago
Location: Fairfax, VA

Re: Question about Simpson's integration.

Post by dirkse »

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
confused student
User
User
Posts: 3
Joined: 5 years ago

Re: Question about Simpson's integration.

Post by confused student »

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

Code: Select all

$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.
Post Reply