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}

Advertisements

Leave a Comment »

No comments yet.

RSS feed for comments on this post.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Blog at WordPress.com.

%d bloggers like this: