# Question:how to use the solution of ODE system for explicit expression and plot it?

## Question:how to use the solution of ODE system for explicit expression and plot it?

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
﻿