Lano2005C.html

> restart;
with(linalg):with(plots):

primitiff0:=proc(pwf,x)
local i,ik,nf,x0,w,pfnai;
ik:=0;
nf:=[];
x0:=0;
pfnai :=0;
for i from 1 to nops(pwf)-2 by 2 do
ik:=eval(pfnai,w=x0);
pfnai := simplify(ik+int(op(pwf)[i+1],x=x0..w));
nf:=[op(nf),subs(s=w,op(pwf)[i]),pfnai];
x0:= rhs(op(pwf)[i]);
#print (op(pwf)[i],pfnai);
od;

ik:=eval(pfnai,w=x0);
pfnai := simplify(ik+int(op(pwf)[i],x=x0..w));
nf:=[op(nf),pfnai];

piecewise(op(subs(w=s,nf)));
end:

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

> N:=30;
L:=1;
M:=1;
g:=1;
rad_spline := 1;

> si := [L*(i-1/2)/(N)$i=1..N];
ui := [phi[i]$i=1..N];
oi := [omega[i]$i=1..N];

> x:=primitiff0(spline(si,map(sin,ui),s,rad_spline),s):
y:=primitiff0(spline(si,map(cos,ui),s,rad_spline),s):
x1:=primitiff0(spline(si,[seq(cos(ui[i])*oi[i],i=1..N)],s,rad_spline),s):
y1:=primitiff0(spline(si,[seq(-sin(ui[i])*oi[i],i=1..N)],s,rad_spline),s):

> #soucet druhych mocnin rychlosti je nad moznosti simplify
vx2vy2 := [op(x)]:
for i in [2*j$j=1..N-2,2*(N-2)+1] do
vx2vy2[i]:=simplify(op(x1)[i]^2+op(y1)[i]^2);
#print (op(x)[i]);
od:
vx2vy2:=piecewise(op(vx2vy2)):
T:=primitiff0(vx2vy2,s):
T:=eval(T,s=L):
T:=simplify(1/2*T);

V:=primitiff0(y,s):
V:=eval(V,s=L):
V:=simplify(-g*M/L*V);

>

> # plne nelinearni Lagrangeovy rovnice
Lt:=subs(map(q->q=q(t),ui),seq(oi[i]=diff(phi[i](t),t),i=1..N),convert(grad(T-V,oi),list)):
ma:=map(simplify,map(diff,Lt,t)):
ma:=subs(seq(diff(phi[i](t),t,t)=epsilon[i],i=1..N),seq(diff(phi[i](t),t)=oi[i],i=1..N),seq(phi[i](t)=phi[i],i=1..N),(ma)):
F:=array(map(combine,simplify(convert(grad(T-V,ui),list)-subs(seq(epsilon[i]=0,i=1..N),ma))));
Q:=matrix(map(combine,jacobian(ma,[seq(epsilon[i],i=1..N)])));
Q:=stackmatrix(Q,F):

>

> #codegen[C](Q,optimized);

>

> fajlik:=fopen("/home/ledvinka/projects/lano/deklarace.c",WRITE):
fprintf(fajlik,"const int StupnuVolnosti = %d;\n",N):
fprintf(fajlik,"const double Phi_0 = %f;\n",1.57):
fprintf(fajlik,"const double dt = %f;\n",1/20):
fprintf(fajlik,"const double tf = %f;\n",20):
fclose(fajlik);

> fajlik:=fopen("/home/ledvinka/projects/lano/rhs.c",WRITE):
for i from 0 to N*6 do
fprintf(fajlik,"double t%d",i*N+1);
for j from 2 to N do
fprintf(fajlik,",t%d",i*N+j);
od;
fprintf(fajlik,";\n");
od:
fclose(fajlik);
codegen[C](Q,optimized,filename="/home/ledvinka/projects/lano/rhs.c");

>

> #READ ALL TIMES
video:=[]:
FN:="/home/ledvinka/projects/lano/vystup.txt";
line:=1:
while (line <> 0) do
line := readline(FN);
if (length(line)>10) then
hodnoty := op(sscanf(line,"%hf"));
nnr:=seq(phi[i]=hodnoty[i+1],i=1..N):
video:=[op(video),plot( [ subs(nnr,x), -subs(nnr,y), s=0..1] ,scaling=constrained,view=[-1..1,-1..0])];
fi;
od:
close(FN);

display(video,insequence=true,scaling=constrained);

> display(video,insequence=true,scaling=constrained,axes=none,color=black);

>
plot( [ subs(nnr,x), -subs(nnr,y), s=0..1] ,scaling=constrained,view=[-1..1,-1..0]);

>

> hodnoty;

>

> with(plots):

Warning, the name changecoords has been redefined

> nnr:=unum(2.0):
nnr:=seq(phi[i]=nnr[i],i=1..N):
plot( [ subs(nnr,x), -subs(nnr,y), s=0..1] ,scaling=constrained,view=[-1..1,-1..0]);

> fps:=20:
video:=[]:
for nt from 0 to eval(20*fps) do
nnr:=unum(nt/fps):
nnr:=seq(phi[i]=nnr[i],i=1..N);
video:=[op(video),plot( [ subs(nnr,x), -subs(nnr,y), s=0..1] ,scaling=constrained,view=[-1..1,-1..0])];
od:

Error, (in rdr) too many function evaluations in rkf45  30004

> display(video,insequence=true,scaling=constrained);
#VIDEO(N=2)

>