Double pendulum model

The state vector: ((((Subscript[θ, 1]),(Subscript[θ, 2])))^T)
The unknown forces: ((((Subscript[F, 1]),(Subscript[F, 2])))^T)

params={g->9.81,m1->1,m2->1,l1-> 1,l2->1};    (* pendulum parameters *)

p1[t_]:={l1 Sin[θ1[t]],-l2 Cos[θ1[t]]}                           (* links end points *)
p2[t_]:={l1 Sin[θ1[t]]+l2 Sin[θ2[t]],-l1 Cos[θ1[t]]-l2 Cos[θ2[t]]}

vg={0,-g};                                                                      (* gravity forces *)
fg1=m1 vg; fg2=m2 vg; 

ft1[t_]:=-F1[t]p1[t]/l1+F2[t](p2[t]-p1[t])/l2   (* tensions *)
ft2[t_]:=-F2[t](p2[t]-p1[t])/l2

eqns={m1 p1''[t]==ft1[t]+fg1,                        (* motion equations *)
m2 p2''[t]==ft2[t]+fg2};