Numerically solving differential equations for multiple initial conditions with Maple

Question

I need to draw a dozen of field lines using Maple, by solving three differential equations. Here's a fully working code for Maple, to be modified :

restart:
with(linalg):
with(DEtools):

Position := [x(t), y(t), z(t)]:
r := sqrt(x(t)^2 + y(t)^2 + z(t)^2):

Mu := [0, 0, 1]:
DipolarField := 3*innerprod(Mu, Position)*Position - Mu*r^2:

Bx := innerprod([1, 0, 0], DipolarField):
By := innerprod([0, 1, 0], DipolarField):
Bz := innerprod([0, 0, 1], DipolarField):

s1 := -50:
s2 := 50:

eq1 := (D(x))(t) = Bx:
eq2 := (D(y))(t) = By:
eq3 := (D(z))(t) = Bz:

x0 := 1.0:
y0 := 0.0:
z0 := 0.0:

InitConditions := x(0) = x0, y(0) = y0, z(0) = z0:
FieldLines := dsolve({eq1, eq2, eq3, InitConditions}, {x(t), y(t), z(t)}, numeric, range = s1..s2, maxfun=0):

plots[odeplot](FieldLines, [x(t), y(t), z(t)], s1..s2, numpoints = 1000, axes = boxed, scaling = constrained, thickness = [3], color = ["Green"]);

This code draws a nice lonely field line. But then I need that code to simultaneously draw 12 field lines, all starting on the equatorial circle, and show them all in the plot. I thus need to change the initial conditions x0 and y0 to something like cos(theta) and sin(theta), and tell Maple to do the numerical integration for theta = 0, Pi/12, 2Pi/12, 3Pi/12, 4Pi/12, ..., 11Pi/12, but I don't know how.

I'm not an experienced user of Maple (I usually use Mathematica), so I dont know how to do this with Maple. Any suggestion ?

Answer 1

You should just have to stick it in a MAPLE do loop, something like:

restart:
with(linalg):
with(DEtools):

Position := [x(t), y(t), z(t)]:
r := sqrt(x(t)^2 + y(t)^2 + z(t)^2):

Mu := [0, 0, 1]:
DipolarField := 3*innerprod(Mu, Position)*Position - Mu*r^2:

Bx := innerprod([1, 0, 0], DipolarField):
By := innerprod([0, 1, 0], DipolarField):
Bz := innerprod([0, 0, 1], DipolarField):

s1 := -50:
s2 := 50:

eq1 := (D(x))(t) = Bx:
eq2 := (D(y))(t) = By:
eq3 := (D(z))(t) = Bz:

theta_list:= [seq(i*Pi/12, i=0..11)]: # YOUR THETA VALUES

COLOUR_LIST:= [blue, black, red, green, yellow, cyan, gold, gray, orange, brown, pink, magenta]: 
# LIST OF COLOURS - CAN BE ANYTHING AS LONG AS YOU HAVE 12 DIFFERENT ONES!

for j from 1 to nops(theta_list) do # THE MAGIC 'DO' LOOP

    theta_val:= theta_list[j]: # VALUE OF THETA FOR EACH j IN LOOP

    COL:= COLOUR_LIST[j]: # ASSIGNS A DIFFERENT COLOUR FOR EACH j IN LOOP

    x0 := cos(theta_val):
    y0 := sin(theta_val):
    z0 := 0.0:

    InitConditions := x(0) = x0, y(0) = y0, z(0) = z0:
    FieldLines := dsolve({eq1, eq2, eq3, InitConditions}, {x(t), y(t), z(t)}, numeric, range = s1..s2, maxfun=0):

    PLT[j]:= plots[odeplot](FieldLines, [x(t), y(t), z(t)], s1..s2, numpoints = 1000, axes = boxed, scaling = constrained, thickness = [3], color = COL); 
    # ASSIGNS A NAME FOR EACH PLOT FOR EACH THETA VALUE 
    # NOTE EACH PLOT IS A DIFFERENT COLOUR DUE TO "color=COL"

od: # END OF j LOOP

plots[display]( seq(PLT[k], k=1..nops(theta_list)) ); 
# DISPLAYS ALL PLOTS FOR ALL THETA VALUES TOGETHER, EACH A DIFFERENT COLOUR

I have tried to annotate the changes. I also added different colours for each plot so that hopefully you can distinguish each curve. If you don't want that, just delete this list and change the color back to color=["Green"] in odeplot. I haven't debugged it as I don't have MAPLE on the machine I am currently using, but hope this helps!