$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;
Any advice on how to proceed is welcome.
Thanks for reading.