# Rcjp's Weblog

## November 7, 2006

### Series Solution Method in maxima

Filed under: physics — rcjp @ 5:56 pm

Just a small script to illustrate series solution

 /*************************************************************************** Solve the Hamiltonian 2 2 2 m q w p H(p, q) = ------- + --- 2 2 m by splitting to two coupled equations d p -- (q(t)) = - dt m d 2 -- (p(t)) = - m q w dt and using series method ****************************************************************************/ kill(p,q,t,m,w)$eq1: diff(q(t), t)= p(t)/m; eq2: diff(p(t), t)= -m*w^2*q(t); assume(w>0)$ atvalue(p(t), t=0, 0)$atvalue(q(t), t=0, 1)$ ans: desolve([eq1, eq2], [q(t), p(t)]); /* And now do it just using power series */ kill(a,b,p,q,t,w,m)$p:taylor(sum(a[i]*t^i,i,0,6),t,0,6)$ q:taylor(sum(b[i]*t^i,i,0,6),t,0,6)$/* Note rearrange eqn=0 */ eq1:diff(q,t)-p/m$ eq2:diff(p,t)+m*w^2*q$coeffs1:makelist(ratcoef(eq1,t,i),i,0,6)$ coeffs2:makelist(ratcoef(eq2,t,i),i,0,6)$eqns:append(cons(subst(0,t,p) = 0,coeffs1),cons(subst(0,t,q) = 1,coeffs2))$ ab:append(makelist(a[i],i,0,6),makelist(b[i],i,0,6))$ans:solve(eqns,ab)$ trunc(subst(ans[1],p)); trunc(subst(ans[1],q)); 

or inside an emaxima session:

 % -*-EMaxima-*- \documentclass{article} \usepackage{a4wide} \usepackage{emaxima} \begin{document} \section*{Series Solution of PDE} $$H(p,q) = {p^2\over 2 m} + {m\over 2} w^2 q^2$$ $$\frac{dq}{dt}=\frac{\partial H}{\partial p} = \frac{p}{m} \qquad{\rm and}\qquad \frac{dp}{dt}=-\frac{\partial H}{\partial q} = -m\omega^2q$$ \begin{maximasession}   kill(p,q,t,m,w)$eq1: diff(q(t), t)= p(t)/m; eq2: diff(p(t), t)= -m*w^2*q(t); assume(w>0)$   atvalue(p(t), t=0, 0)$atvalue(q(t), t=0, 1)$      ans: desolve([eq1, eq2], [q(t), p(t)]); \maximaoutput* \end{maximasession} To do that with series expansions only: \begin{maximasession}   p: taylor(sum(a[i]*t^i,i,0,inf),t,0,6)$q: taylor(sum(b[i]*t^i,i,0,inf),t,0,6)$    eq1: diff(q, t) - p/m$eq2: diff(p, t) + m*w^2*q$   coeffs1: makelist (ratcoeff(eq1, t, i), i , 0 , 6)$coeffs2: makelist (ratcoeff(eq2, t, i), i , 0 , 6)$   /* add in boundary conditions p=0 and q=1 at t=0 */   eqns: append(cons(subst(0,t,p)=0, coeffs1), cons(subst(0,t,q)=1, coeffs2))$coeffs:append(makelist(a[i],i,0,6), makelist(b[i],i,0,6))$   ans: solve(eqns, coeffs)\$   trunc(subst(ans[1], p));   trunc(subst(ans[1], q)); \maximaoutput* \end{maximasession} \end{document}