Fancy2.mws

**plot embellishments**

**more general plotting of dsolve results (Perry's magic)**

**automated solving of N equations in N unknowns**

**charge distribution over a wire**

**rotating mass-spring**

**restart; with(plots): **

The contour plot is a special case of plot3d, which shows contours. Any plot3d graph can be manipulated (oriented) to give contour plots.

**f:=1/sqrt((x-a)^2+(y-b)^2);
contourplot(subs(a=-1,b=1,f)-subs(a=1,b=2,f),x=-2..1,y=0..3,contours=20,grid=[60,60],view=-5..5)**;

The 'view' instruction limits the range of the z-direction, and the 'grid' command gives the grid on which the contour is based.

The 3d view is a manipulated view of the contour plot.

Here is a problem where a 40x40 is set up via a loop and then solved
and the results graphed. **This problem is one model of
how charge accumulates at sharp points**, in this case, at the
ends of a wire.

This calculates the approximate surface charge on a long thin wire, where the wire radius is small compared to its overall length.

It determines the charges on each subsection of the wire (solving N equations in N unknowns), by specifying that the potential on the axis of each hollow cylindrically-shaped element be the same V.

**restart; with(plots):**

Potential due to a charged section of a cylinder when the ends of the cylinder are located at Z1 and Z2

**V:=Int(Int(k*rho*2*Pi*r/sqrt(r^2+z^2),r=0..R),z=Z1..Z2);
V:=value(V);** # Evaluate the double integral.

**L:=5;** #
The length of a section

**Vs:=1;** # The potential measured at
each section

Find the potential at different segment distances away from a charged segment.

**d:=(dist)->evalf(subs({Z1=-L/2+L*dist,Z2=L/2+L*dist},Vsub));**

Make a function for the potential of the section at "dist" sections away from the charged section with rho="charge".

**Vsec:=(dist,charge)->subs(rho=charge,d(dist));**

Set up 4 equations to find the 4 unknown charges on the 4 sections.

**V1:=Vs=Vsec(0,rho1)+Vsec(1,rho2)+Vsec(2,rho3)+Vsec(3,rho4);
V2:=Vs=Vsec(1,rho1)+Vsec(0,rho2)+Vsec(1,rho3)+Vsec(2,rho4); V3:=Vs=Vsec(2,rho1)+Vsec(1,rho2)+Vsec(0,rho3)+Vsec(1,rho4);
V4:=Vs=Vsec(3,rho1)+Vsec(2,rho2)+Vsec(1,rho3)+Vsec(0,rho4);**

**fsolve({V.(1..4)},{rho.(1..4)});**
# Solve four equations for four unknowns

Now do the same problem with more segments in the cylinder.

**NumSeg:=40;** # Number of segments
in cylinder.

Use a loop to create the equations.

**for m from 1 to NumSeg do
V.m:=Vs=convert([seq(Vsec(abs(m-n),rho.n),n=1..NumSeg)],`+`): od:**

This loop creates NumSeg equations. In each equation
Vs is equated to a sum of items.

This sum is converted (via '+' ) from a sequence of items, created from
the Vsec function..

**sol:=fsolve({V.(1..NumSeg)},{rho.(1..NumSeg)});**
# Solve the 40 equations for the 40 unknowns.

**assign(sol); **# Assign the solutions
to p1-p40

**plot([seq([n,rho.n],n=1..NumSeg)],title=`Charge
distribution across the bar`); **# plot the results.

Here is complex problem solved via Lagrange's equations using a numeric version of dsolve. Then some creative functional assignments are made which allow the results to be plotted.

Two equal masses are connected by a spring, and are rotating about their CM.

rotating mass-spring system in Lagrangian formalism, patterned after proj.ms Maple seems to need a lot of shifting around in the variables.

**restart;with(plots):
alias(r=r(t),phi=phi(t)):
x1:=r/2*cos(phi): y1:=r/2*sin(phi): x2:=-x1:y2:=-y1: KE:=simplify(m/2*(diff(x1,t)^2+diff(y1,t)^2+diff(x2,t)^2+diff(y2,t)^2));
PE:=k/2*(r-r0)^2;
L:=KE-PE;** # lagrangian

**valuesin:={diff(r,t)=rdot,diff(phi,t)=phidot,r=ra,phi=phia};
Lsvel:=subs(valuesin,L);** # L with dummy variables

Generalized momenta w/ symbols (pj = dLdqjdot)

**pr:=diff(Lsvel,rdot); pphi:=diff(Lsvel,phidot); **

**dLdr:=diff(Lsvel,ra);dLdphi:=diff(Lsvel,phia); **# dL/dqj

**valuesback:={rdot=diff(r,t),phidot=diff(phi,t),ra=r,phia=phi};**

**pra:= subs(valuesback,pr);
pphia:=subs(valuesback,pphi); **# generalized p's 'live'

dLdphia:=subs(valuesback,dLdphi);

**eqr:= diff(pra,t)-dLdra = 0; **# lagrange's eqn in r.

**eqphi:=diff(pphia,t)-dLdphia = 0;** # ditto for phi.

**eqra:= subs(r=a(t),phi=b(t),m=1,r0=1,k=2,eqr);
eqphia:=subs(r=a(t),phi=b(t),m=1,r0=1,k=2,eqphi);** # new set of named
vars

The square brackets are important in the following dsolve command to keep the solutions in correct order

**sol:=dsolve({eqra,a(0)=3/2,D(a)(0)=0,eqphia,b(0)=0,D(b)(0)=10},[a(t),b(t)],numeric);
sol(1); **# take a look at the solution when t=1 sec

**rp:= z -> rhs(sol(z)[2]); phip:= z -> rhs(sol(z)[4]);**

magic assignments (Perry worked this out)

**h:=plot(rp,0..3,color=blue):
j:=plot(phip,0..3,color=red):
display({h,j},title=` r(blue) and phi(red) vs t`);**

**plot([rp,phip,0..12],coords=polar,title=`polar plot of motion`);**

End of Fancy2.mws