i got the solution of a ODE system (dsn), the results are time dependent variable theta[i](t),and its deriavtives diff(theta[i](t),t). and i have one set of expressions, which are explicit form of theta[i](t),and diff(theta[i](t),t), such as ss[1]:=fx[1]=sin(theta[1](t))*3+cos(theta[1](t))*diff(theta[2](t),t). I want to put the solution i got into the explicit expression of fx[i], and then polt them, how can I do it? i use eval(ss[1],dsn)but it does not work. here is the code:
> restart:
interface(warnlevel=0):
with(SolveTools):with( plots ): with( plottools ):
with( DEtools ): with( PDEtools ):
g:=9.81:
x[0][2]:=0:
y[0][2]:=0:
#initialization function of each beam,in this function, we should give 3 variables, n=the number of this #beam, la=the length of this beam, ma=the mass of this beam
beams:=proc(n,la,ma)
global l,m,x,y:
l[n]:=la: #length of beam n
m[n]:=ma: #mass of beam n
#initialize the positions of each beam
x[n][2]:=0: #the x position of the end point of each beam
y[n][2]:=0: #the x position of the end point of each beam
x[n][1]:=0: #the y position of the center point of each beam
y[n][1]:=0: #the y position of the center point of each beam
end proc;
#define the forces that applied on the end of beam n
Force:=proc(n,fm,fn)
global fx,fy:
fx[n]:=fm:
fy[n]:=fn:
end proc:
> beams(1,1,1):
beams(2,2,1):
beams(3,3,1):
beams(4,4,1);
> Force(2,80,100);
> la:=proc(nn)
global x,y,q,vy,vx,T,lg,lm,lp,L,Ll,Lk,ax,ay:
n:=nn: #transfer the function varibale to local function
# generate the position of the beams, nn is the number of beams
# x,y is the position vector, q is the angle vector
for i from 1 to n do:
x[i][2]:=l[i]*cos(q[i](t))+x[(i-1)][2]: #x[i][2],y[i][2] is the positon of the end of beam i
y[i][2]:=l[i]*sin(q[i](t))+y[(i-1)][2]:
x[i][1]:=x[i][2]-0.5*l[i]*cos(q[i](t)): #x[i][1],y[i][1] is the center of beam i
y[i][1]:=y[i][2]-0.5*l[i]*sin(q[i](t)):
end do:
#generate the velocity of the center of the beams
for i from 1 to n do:
vx[i]:=diff(x[i][1],t):
vy[i]:=diff(y[i][1],t):
ax[i]:=diff(x[i][1],t$2):
ay[i]:=diff(y[i][1],t$2):
end do:
#generate the kinetic energy of the system
T:=0:
for i from 1 to n do:
T:=T+0.5*m[i]*(vx[i]^2+vy[i]^2)+1/24*m[i]*l[i]^2*diff(q[i](t),t)^2:
end do:
# substitue the angle and the angle velocity with time independent variables, for later diffrenciation
for i from 1 to 4 do:
T:=subs(diff(q[i](t),t)=omega[i],T); # angle velocity(q(t) is angle function)
T:=subs(q[i](t)=theta[i],T): #angle (q(t) is angle function))
end do;
#generate the lagrange equations
#lg,lm,lp,L,Ll,Lk,those variables is the intermedia variables for Lagrange equation, explained detailly below
for i from 1 to n do:
lg[i]:=diff(T,omega[i]): #derivation of "T" with respect to angle velocity
#substitute the angle and angle velocity in each lg[i] with time function, for later derivation over time
for j from 1 to n do:
lg[i]:=subs({theta[j]=theta[j](t),omega[j]=omega[j](t)},lg[i]):
end do;
end do:
#---------------------------------------------------------
for p from 1 to n do:
lp[p]:=diff(lg[p],t) :#derivation of "lg" term over time,make lp the first term in lagrange eq.
end do:
#--------------------------------------------------------
for o from 1 to n do:
lm[o]:=diff(T,theta[o]): #derivation of "T" over angle, which is the 2nd term in the eq.
for k from 1 to n do:
#substitute the angle and angle velocity in each lm[o] with time function, for later derivation over time
lm[o]:=subs({theta[k]=theta[k](t),omega[k]=omega[k](t)},lm[o]):
end do:
end do:
#---------------------------------------------------------
for k from 1 to n do:
#L is the lagrange eq. with the virtual displacement of angle
#Ll is the lagrange eq. with the virtual displacement of x
#Lk is the lagrange eq. with the virtual displacement of y
#unlike L,Ll and Lk are written directly,no deduction by programming
L[k]:=-lp[k]+lm[k]+(fy[k]+0.5*m[k]*g)*cos(theta[k](t))*l[k]-fx[k]*sin(theta[k](t))*l[k]=0:
#L[k]:=-lp[k]+lm[k]+(fy[k]+0.5*m[k]*g)*cos(theta[k](t))*l[k]-fx[k]*sin(theta[k](t))*l[k]-(fy[k]+0.5*m[k]*g)*theta[k](t)*sin(theta[k]##(t))*l[k]-fx[k]*cos(theta[k](t))*l[k]*theta[k](t)=0:
#because diff(T,x)=Diff(T,omega)*Diff(omega,x)
#Ll[k]:=lp[k]*(1/l[k])*cos(theta[k](t))*omega[k](t)/sin(theta[k](t))^2-lm[k]/(l[k]*sin(theta[k](t)))-fx[k]:
#because diff(T,y)=Diff(T,omega)*Diff(omega,y)
#Lk[k]:=lp[k]*(1/l[k])*sin(theta[k](t))*omega[k](t)/cos(theta[k](t))^2+lm[k]/(l[k]*cos(theta[k](t)))-fy[k]-m[#k]*g:
Ll[k]:=fx[k-1]-fx[k]=m[k]*ax[k]:
Lk[k]:=fy[k-1]-fy[k]-0.5*m[k]*g=m[k]*ay[k]:
end do:
for s from 1 to n do:
for r from 1 to n do:
L[s]:=subs(diff(omega[r](t),t)=diff(theta[r](t),t,t),L[s]):
#Lk[s]:=subs(diff(omega[r](t),t)=diff(theta[r](t),t,t),Lk[s]):
Lk[s]:=subs(q[r](t)=theta[r](t),Lk[s]):
Ll[s]:=subs(q[r](t)=theta[r](t),Ll[s]):
#Ll[s]:=subs(diff(omega[r](t),t)=diff(theta[r](t),t,t),Ll[s]):
L[s]:=subs(omega[r](t)=diff(theta[r](t),t),L[s]):
#Lk[s]:=subs(omega[r](t)=diff(theta[r](t),t),Lk[s]):
#Ll[s]:=subs(omega[r](t)=diff(theta[r](t),t),Ll[s]):
L[s]:=combine(L[s],trig):
Ll[s]:=combine(Ll[s],trig):
Lk[s]:=combine(Lk[s],trig):
end do:
end do:
end proc;
> la(2):
solve_la:=proc(nn,inc)
global L,Lk,Ll,dsys:dsys2,ini,dsn:
ini:=inc:
n:=nn:
dsys:={seq(Lk[i],i=1..n),seq(Ll[i],i=1..n)}:
ss:=solve(dsys,{seq(fx[i],i=0..n-1),seq(fy[i],i=0..n-1)}):
dsys2:=eval({seq(L[i],i=1..n)},ss);
dsn := dsolve({op(dsys2)} union ini,numeric,implicit=true);
end proc;
n:=2:
j:=n-1:
#Lk[1];
#Lk[2];
#Ll[1];
#Ll[2];
#L[1];
#L[2];
#dsys:={Lk[2],Lk[1],Ll[1],Ll[2],L[1],L[2]};
dsys:={seq(Lk[i],i=1..n),seq(Ll[i],i=1..n)};
#dsys:={Lk[2],Lk[1],Ll[1],Ll[2]};
> ss:=solve(dsys,{seq(fx[i],i=0..n-1),seq(fy[i],i=0..n-1)});
> dsys2:=eval({seq(L[i],i=1..n)},ss);
> ini := {theta[1](0) = 1/6*Pi, theta[2](0) = 1/3*Pi, D(theta[1])(0) = 0, D(theta[2])(0) = 0};
dsn := dsolve({op(dsys2)} union ini,type=numeric,output=listprocedure);
#here is the problem, i want to put dsn into #ss, then get the explicit expression of #fx,and fy, then plot them, how?
eval(ss,dsn);
Error, invalid input: cos expects its 1st argument, x, to be of type algebraic, but received proc (t)::; local res, data, solnproc, outpoint, `theta[2](t)`; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 14; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](t) else outpoint := evalf(t) end if; data := module () local Data; export Get, Set; end module; solnproc := data:-data("soln_procedure"); if not type(outpoint, 'numeric') then if member(t, ["start", 'start', "method", 'method', "left", "right", "leftdata", "rightdata", "enginedata"]) then if member(t, ['start', 'method']) then res := solnproc(convert(t, 'string')) else res := solnproc(t) end if; if type(res, 'array') then return eval(res, 1) elif res <> "procname" then return res end if elif member(t, ["last", "initial"]) then res := solnproc(t); if type(res, 'list') then return res[4] end if elif member(t, ['last', 'initial']) then res := solnproc(convert(t, 'string')); if type(res, 'list') then return res[4] end if elif type(outpoint, `=`) and member(lhs(outpoint), ["initial", 'initial']) then if type(rhs(outpoint), 'list') then res := solnproc("initial" = [0, op(rhs(outpoint))]) else res := solnproc("initial" = [4, rhs(outpoint)]) end if; if type(res, 'list') then return res[4] end if elif outpoint = "solnprocedure" then return eval(solnproc) elif outpoint = "sysvars" then return data:-data("sysvars") end if; if procname <> unknown then return 'procname'(t) else `theta[2](t)` := pointto(data:-data("soln_procedures")[4]); return '`theta[2](t)`'(t) end if end if; try res := solnproc(outpoint); res[4] catch: error end try end proc
> plots[odeplot](dsn,[[t,theta[1](t)],[t,theta[2](t)]],0..1);
1.6+ BBBB
+ B B BB
1.4+ AA B BB A BB BBB AA
+ AAA AA B B AA AAA BB BB AA AA
1.2+ AA AAB BB AA AA B BBA AA
*B AA BA B AA AA B *B AA BBBBBB
1+ BB AA BB A B*A AA B A BB A*B BB
+ B*A B A A* A* AA BB BAA BB
+ A*B B A ABB B*A A BBBBB AA **
0.6+ A BB BB AA AA B B AA A AA AA
*AA BBBB AAAA BB BB AAAA AAAA
0.4+ B B
+ BB BB
0.2+ BB BB
-+--+--+--+--+--+-+--+--+--+--+--+****+--+--+--+--+--+--+-+--+--+--+--+--+-
0 0.2 0.4 0.6 0.8 1
AAAAAAAAAAAAAAAA theta[1]
BBBBBBBBBBBBBBBB theta[2]
>
This post generated using the online HTML conversion tool
Download the original worksheet
