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.
Vsub:=subs({k=1,R=1},V); # Set the value of the k constant to 1, and the radius of the cylinder to 1.

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'
dLdra:=subs(valuesback,dLdr);
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