Dr. Patrick T

## 2083 Reputation

13 years, 301 days

## DEtools[DEplot]...

Hi Spring,

You may find useful hints under DEtools[DEplot]

The difficulty with your problem is that you have 4D system but you would like a 2D phase plot. There exist built-in routines to create 2D phase plots for 2D systems, namely dfieldplot and phaseportrait.

A fun example of 2D phase plot is given in the Help menu:

```vdP := [diff(x(t),t)=10*(y(t)-x(t)^3/3+x(t)),diff(y(t),t)=-1/10*x(t)];
DEtools[DEplot](vdP,[x(t),y(t)],t=0..20,x=-3..3,y=-2..2,[[x(0)=0.2,y(0)=0.2]],
title=`van der Pol oscillator`,dirfield=400,arrows=comet,
size=magnitude,numpoints=505,linecolor=blue,animatecurves=true);

```

There is something for 3D systems and 3D phase plots, namely fieldplot3d, but I have found it difficult to make nice plots with it. You may have better luck.

At any rate you have a 4D system so none of these methods is directly applicable.

One second-best approach is to look for a 2D approximation of your system near some point of interest.

Another second-best is to use a color scheme to show the direction of motion. It doesn't have arrows though.

```DEtools[phaseportrait]([D(x)(t)=y(t)-z(t),D(y)(t)=z(t)-x(t),D(z)(t)=x(t)-y(t)*2],
[x(t),y(t),z(t)],t=-2..2,[[x(0)=1,y(0)=0,z(0)=2]],stepsize=.05,
scene=[z(t),x(t)],linecolour=sin(t*Pi/2),method=classical[foreuler]);
```

As jakubi has suggested you may have to hard-code your own arrows based on the direction field in the plane of interest. If you manage to do something, I'd love to hear about it, so please do post your solution here once you have it.

Quite often, the more code you post, the more useful feedback you will get.

## fantastic...

A great thank you Robert. It's EXACTLY what I was after. It runs great on my machine and will save a lot of copy-pasting. It will also help me have consistent options on all my plots throughout the worksheet.

A great thank you Joe. I'll remember this "try/finally" trick. Very useful.

And by the way those were two extremely quick answers! I haven't even had to interrupt my work waiting for the solution. A godsend.

## maple worksheet problem...

This is the header of your worksheet. It's full of html formatting. Can you post a pure Maple version?

<!DOCTYPE html
PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">

## Maple 11.02 supports tickmarks,thickness...

A related thread for future reference.

As jakubi has shown in the link below, the following works in Maple 11.02:

plot(cos(x), x=-2*Pi..2*Pi, axis=[tickmarks=[10, color=blue, thickness=3]]);

See:

http://www.mapleprimes.com/forum/tickmarkshowsetspacingandthickness#comment-26122

## Maple 11 rules!...

that's a beautiful graph, now I want blue tickmarks more than ever!

and here's another case where you wish you could type setversion(11) or something like that into Maple 13.

Thanks David Clayworth, Thanks jakubi, I'll try to get ahold of V.11 (I have the disks somewhere, but where...)

## language...

the meaning of "without approximations" seems to be improper here: numerical solutions of ode systems (such as Runge-Kutta schemes) are approximations. Did the prof actually write that or is it your recollection?

you can copy-paste the Maple code into the window. If it's long, try to break the execution blocks and to insert linebreaks to make it readable.

Before copy-pasting, remove the output  (Edit->Remove Output) to remove images and the likes.

posting a worksheet is also very easy, once you've posted it you can just copy-paste the url link to the worksheet. Go into "my files" on the lhs menu.

## nonlinear...

It's highly unlikely that there exist transformations that can make this ode linear. Nonlinear it is, nonlinear it will certainly remain. You just don't have a choice, you have to solve it numerically, assuming a solution does exist.

>That nonlinear DE in a(r) can be solved by using 'equilibrium iteration method or equilibrium Newton's method' numerically.

why don't you show the Maple code you've used to do that?

## taylor approximate polylog...

If you're so confident that there's a solution, then why don't you start with an approximation of polylog in terms of polynomials. Check the definition of the polylog (in Maple's Help menu), take a few terms of the expansion and see if that gives anything.

Also note that if 3/2 is replaced by something else (say by 1), the polylog function sometimes simplifies to better known functions.

But if I were you I would listen very, very carefully to Robert's advice.

## to get you started...

Help:

dsolve[numeric,bvp],

you've got the following problems: some undefined constants, the polylog function, a bvp defined by a second-order ode.

to get started, I tried this:

> C:=1: C2:=1: C3:=1:
> ode := diff(y(r), r, r) +2*(diff(y(r), r))/r + C*polylog(3/2, -exp(C2*(C3-y(r))));
> R := 100:
> bcs := D(y)(0)=0, y(R)=0:
> dsolve({ode, bcs}, {y(r)}, range=0..R, type=numeric):
Error, (in dsolve/numeric/bvp) system is singular at left endpoint, use midpoint method instead

don't know how to proceed from here, sorry.

you could also transform the bvp into a set of 2 odes, but I don't know if that would help.

## gridlines works as expected...

works:

```plot(cos(x), x=-2*Pi..2*Pi,
axis = [gridlines = [10, color = blue, thickness = 3]]);

```

doesn't work:

```plot(cos(x), x=-2*Pi..2*Pi,
axis = [tickmarks = [10, color = blue, thickness = 3]]);

```

## cool colors...

by combining:

style = patchnogrid,
lightmodel = light1,

EDIT:

Just noticed: the cool effect is in Maple Classic. In Maple Standard, it looks totally different...

## another option...

I think jakubi's suggestion can be obtained with the additional option:

## plot3d options : lightmodel, ambientligh...

Check this:

plot3d options : lightmodel, ambientlight, color

Perhaps lightmodel=light1 does what you're looking for. Play around with this and let us know which options you found best.

Notice how line breaks (SHIFT+RETURN) make it easier to read your code.

And creating a separate "execution group" for the plot allows you to play with your options quicker, in that you don't have to recompute everything each time.

I have inserted extra options that are available by right-clicking, in the way suggested by jakubi above, and which you can also manually force. See if you can set the color to Z (Hue) (I couldn't off the top of my head just now).

> restart:
with(LinearAlgebra):
with(plots):

> N := 10:
n1 := 1.5:
n2 := 3.4:
n3 := 3.5:
d1 := 0.16e-6:
d2 := d1:
d3 := d1:
c := 0.3e9:
theta2 := arcsin(n1*sin(theta1)/n2):
theta3 := arcsin(n2*sin(theta2)/n3):
`#mover(mi("n1"),mi("~"))` := n1*cos(theta1):
`#mover(mi("n2"),mi("~"))` := n2*cos(theta2):
`#mover(mi("n3"),mi("~"))` := n3*cos(theta3):
Lambda := d1+d2+d3:
`#mover(mi("n"),mo("&uminus0;"))` := (n1*d1+n2*d2+n3*d3)/Lambda:
`?b` := c/(`#mover(mi("n"),mo("&uminus0;"))`*(2*Lambda)):
`evalf/constant/?B` := proc () options operator, arrow:
`?b` end proc:
phi1 := 2*Pi*n1*d1*nu*cos(theta1)/c:
phi2 := 2*Pi*n1*d1*nu*cos(theta2)/c:
phi3 := 2*Pi*n1*d1*nu*cos(theta3)/c:
M111 := (`#mover(mi("n2"),mi("~"))`+`#mover(mi("n1"),mi("~"))`)*exp(-I*phi1):
M112 := (`#mover(mi("n2"),mi("~"))`-`#mover(mi("n1"),mi("~"))`)*exp(I*phi1):
M121 := (`#mover(mi("n2"),mi("~"))`-`#mover(mi("n1"),mi("~"))`)*exp(-I*phi1):
M122 := (`#mover(mi("n2"),mi("~"))`+`#mover(mi("n1"),mi("~"))`)*exp(I*phi1):
M211 := (`#mover(mi("n3"),mi("~"))`+`#mover(mi("n2"),mi("~"))`)*exp(-I*phi2):
M212 := (`#mover(mi("n3"),mi("~"))`-`#mover(mi("n2"),mi("~"))`)*exp(I*phi2):
M221 := (`#mover(mi("n3"),mi("~"))`-`#mover(mi("n2"),mi("~"))`)*exp(-I*phi2):
M222 := (`#mover(mi("n3"),mi("~"))`+`#mover(mi("n2"),mi("~"))`)*exp(I*phi2):
M311 := (`#mover(mi("n1"),mi("~"))`+n3)*exp(-I*phi3):
M312 := (`#mover(mi("n1"),mi("~"))`-n3)*exp(I*phi3):
M321 := (`#mover(mi("n1"),mi("~"))`-n3)*exp(-I*phi3):
M322 := (`#mover(mi("n1"),mi("~"))`+n3)*exp(I*phi3):
M1 := (Matrix(2, 2, {(1, 1) = M111, (1, 2) = M112, (2, 1) = M121, (2, 2) = M122}))/(2*`#mover(mi("n2"),mi("~"))`):
M2 := (Matrix(2, 2, {(1, 1) = M211, (1, 2) = M212, (2, 1) = M221, (2, 2) = M222}))/(2*`#mover(mi("n3"),mi("~"))`):
M3 := (Matrix(2, 2, {(1, 1) = M311, (1, 2) = M312, (2, 1) = M321, (2, 2) = M322}))/(2*`#mover(mi("n1"),mi("~"))`):
Mcell := `.`(`.`(M3, M2), M1):
t := 1/Mcell[2, 2]:
r := t*Mcell[1, 2]:
h := abs(Re(Mcell[2, 2])):
R := abs(r)^2:
`&varphi;` := arccos(Re(Mcell[2, 2])):
`&varphi;i` := arccosh(abs(Re(Mcell[2, 2]))):
`?N1` := sin(N*`&varphi;`)/sin(`&varphi;`):
`?N2` := sinh(N*`&varphi;i`)/sinh(`&varphi;i`):
RN1 := `?N1`^2*R/(1-R+`?N1`^2*R):
RN2 := `?N2`^2*R/(1-R+`?N2`^2*R):
f := piecewise(h <= 1, RN1, RN2):

> implicitplot3d(z-f, nu = 0 .. 6*`?b`, theta1 = 0 .. (1/2)*Pi, z = 0 .. 1,
tickmarks = [default, [seq((1/18)*i*Pi = (10*i)^`&ofr;`, i = 1 .. 9)], default],
resolution = 10000,
axes = boxed,
style = patch,
orientation = [30,45],
lightmodel = light1,
labels = ["?", "?", "R"]);

## tickmarks suboptions documentation...

thanks jakubi,

the tickness was really the option I cared about (color was added to test multiple options).

Any Maple developers out there can confirm whether there's a problem with the documentation?

thanks!

plot(cos(x), x=-2*Pi..2*Pi, axis=[tickmarks=[10, thickness=3]]);

## good work woodwise...

good work! I've long dreamed of being able to toggle styles, but gave up some time ago ... sob.

﻿