## Polarplot warning for a semistable limit cycle...

Hey, i'm trying do demonstrate that a nonlinear system has a semistable limit cycle but i get a warning at the plot command saying "Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct" and i dont understand it. So i wonder if someone here could help me?

restart; with(PDEtools); with(plots);
eq1 := diff(x(t), t) = x(t)*(x(t)^2+y(t)^2-1)^2-y(t);
2
d              /    2       2    \
--- x(t) = x(t) \x(t)  + y(t)  - 1/  - y(t)
dt
eq2 := diff(y(t), t) = y(t)*(x(t)^2+y(t)^2-1)^2+x(t);
2
d              /    2       2    \
--- y(t) = y(t) \x(t)  + y(t)  - 1/  + x(t)
dt
tr := {x(t) = r(t)*cos(theta(t)), y(t) = r(t)*sin(theta(t))};
{x(t) = r(t) cos(theta(t)), y(t) = r(t) sin(theta(t))}
eq1b := dchange(tr, x(t)*eq1+y(t)*eq2, [r(t), theta(t)], simplify);
/ d      \       2 /        4         2\
r(t) |--- r(t)| = r(t)  \1 + r(t)  - 2 r(t) /
\ dt     /
eq1b := expand(eq1b/r(t));
d                    5         3
--- r(t) = r(t) + r(t)  - 2 r(t)
dt
eq2b := dchange(tr, y(t)*eq1-x(t)*eq2, [r(t), theta(t)], simplify);
2 / d          \        2
-r(t)  |--- theta(t)| = -r(t)
\ dt         /
eq2b := simplify(eq2b/(-r(t)^2));
d
--- theta(t) = 1
dt
sol1 := dsolve({eq1b, r(0) = r[0]}, r(t));
/      /  /     2  \
|      |  | r[0]   |          2     2
r(t) = exp|RootOf|ln|--------| (exp(_Z))  r[0]
\      \  \r[0] - 1/

2     2
- ln(r[0] + 1) (exp(_Z))  r[0]

/             2\
|(exp(_Z) - 1) |          2     2            2        2
- ln|--------------| (exp(_Z))  r[0]  + (exp(_Z))  _Z r[0]
\ exp(_Z) - 2  /

/     2  \
2     2         | r[0]   |             2
+ 2 (exp(_Z))  r[0]  t - 2 ln|--------| exp(_Z) r[0]
\r[0] - 1/

2
+ 2 ln(r[0] + 1) exp(_Z) r[0]

/             2\
|(exp(_Z) - 1) |             2                    2
+ 2 ln|--------------| exp(_Z) r[0]  - 2 exp(_Z) _Z r[0]
\ exp(_Z) - 2  /

/     2  \
2       | r[0]   |          2
- 4 exp(_Z) r[0]  t - ln|--------| (exp(_Z))
\r[0] - 1/

/             2\
2     |(exp(_Z) - 1) |          2
+ ln(r[0] + 1) (exp(_Z))  + ln|--------------| (exp(_Z))
\ exp(_Z) - 2  /

/     2  \
2                   2       | r[0]   |
- (exp(_Z))  _Z - 2 t (exp(_Z))  + 2 ln|--------| exp(_Z)
\r[0] - 1/

/             2\
|(exp(_Z) - 1) |
- 2 ln(r[0] + 1) exp(_Z) - 2 ln|--------------| exp(_Z)
\ exp(_Z) - 2  /

2                                    2
- (exp(_Z))  + 2 _Z exp(_Z) + 4 t exp(_Z) + r[0]  + 2 exp(_Z)

\\
||
- 1|| - 1
//
sol1 := simplify(sol1);
/      /   /     2  \
|      |   | r[0]   |               2
r(t) = exp|RootOf|-ln|--------| exp(2 _Z) r[0]
\      \   \r[0] - 1/

2
+ ln(r[0] + 1) exp(2 _Z) r[0]

/             2\
|(exp(_Z) - 1) |               2                    2
+ ln|--------------| exp(2 _Z) r[0]  - exp(2 _Z) _Z r[0]
\ exp(_Z) - 2  /

/     2  \
2         | r[0]   |             2
- 2 exp(2 _Z) r[0]  t + 2 ln|--------| exp(_Z) r[0]
\r[0] - 1/

2
- 2 ln(r[0] + 1) exp(_Z) r[0]

/             2\
|(exp(_Z) - 1) |             2                    2
- 2 ln|--------------| exp(_Z) r[0]  + 2 exp(_Z) _Z r[0]
\ exp(_Z) - 2  /

/     2  \
2       | r[0]   |
+ 4 exp(_Z) r[0]  t + ln|--------| exp(2 _Z)
\r[0] - 1/

/             2\
|(exp(_Z) - 1) |
- ln(r[0] + 1) exp(2 _Z) - ln|--------------| exp(2 _Z)
\ exp(_Z) - 2  /

/     2  \
| r[0]   |
+ exp(2 _Z) _Z + 2 t exp(2 _Z) - 2 ln|--------| exp(_Z)
\r[0] - 1/

/             2\
|(exp(_Z) - 1) |
+ 2 ln(r[0] + 1) exp(_Z) + 2 ln|--------------| exp(_Z)
\ exp(_Z) - 2  /

2
+ exp(2 _Z) - 2 _Z exp(_Z) - 4 t exp(_Z) - r[0]  - 2 exp(_Z)

\\
||
+ 1|| - 1
//
sol2 := dsolve({eq2b, theta(0) = theta[0]}, theta(t));
theta(t) = t + theta[0]
theta[0] := (1/4)*Pi;
1
- Pi
4
plot1 := polarplot([subs(r[0] = .1, rhs(sol1)), rhs(sol2), t = 0 .. 10], color = red);
Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct
plot2 := polarplot([subs(r[0] = 2, rhs(sol1)), rhs(sol2), t = 0 .. 10], color = blue);
Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct
display({plot1, plot2}, scaling = constrained, tickmarks = [4, 3], view = [-2 .. 2, -2 .. 2]);

## computing the wronskian with Maple gave nothing...

I am tying to compute the wronskian of a fourth order DE: y=C1e2x+ C2e-x +C3xe-x+ C4x2e-x Here's what I did:

with(VectorCalculus):
with(LinearAlgebra):

Determinant(Wronskian([e^(2*x), e^(-x), xe^(-x), x^2*e^(-x)], x)):

which gave nothing.

AJ

## trouble with Levi-Civita symbol from physics-packa...

Hi, there!

I try to simplify the next expression:

Simplify(LeviCivita[4, sigma, lambda, rho]*LeviCivita[4, xi, eta, mu]*g_[rho, mu]*qp[sigma]*q[lambda]*qp[xi]*q[eta])

In reality it is incorrect answer, because indices must run over 1,2,3 but not 1,2,3,4!

SumOverRepeatedIndices(%) confirms that maple mistakes.

my preamble is:

with(Physics)

Setup(mathematicalnotation = true); Coordinates(X);

Setup(spaceindices = lowercaselatin)

Setup(tensors = q[mu](X))

PDEtools:-declare(q(X))

Setup(tensors = qp[mu](X))

PDEtools:-declare(qp(X))

## Emacs Maple Mode...

I use Ubuntu 14.04 and X, not the desktop.  I use emacs/maple 2016.

GNU Emacs 25.1.2 (x86_64-unknown-linux-gnu, X toolkit, Xaw scroll bars)
of 2017-03-1

;;; maplev.el --- Maple mode for GNU Emacs

;; Authors:    Joseph S. Riel <joer@k-online.com>
;;             and Roland Winkler <Roland.Winkler@physik.uni-erlangen.de>
;; Time-stamp: "2003-10-09 22:49:16 joe"
;; Created:    June 1999
;; Version:    2.155
;; Keywords:   Maple, languages
;; X-URL:      http://www.k-online.com/~joer/maplev/maplev.html
;; X-RCS:      $Id: maplev.el,v 1.14 2006-06-02 14:02:38 joe Exp$

I use emacs/maple mode with maple 2016.  Quite often, emacs looses connection with the maple server.  I do  not remember this happening or maybe not as often, with earlier versions of maple.

After using maple/emacs, I started xmaple.  After a few expression evaluations, the maple server stopped.  Restarting xmaple and repeating the expression evaluations many times, I do not get the crash.  So, this appears to be a difficulty with external connections to the maple server.

Does a later version of maple mode exist?

## Plot an undefined name produces a line...

I expected plot with an undefined name to do nothing, but,

plot(asdf);

actually plots y=x!

## Right click menu only some times works...

When I right click on an expression, I only some times get the context menu (under Maple 2016 in document mode, in Win 10).

Any suggestions - Anyone else having this problem? - I know how to use the long forms, but this is for students - for a brief exposure to Maple.

## Error, (in pdsolve/numeric/process_IBCs) initial/b...

hi

how i can remove this error?

thanks

Error, (in pdsolve/numeric/process_IBCs) initial/boundary conditions can only contain derivatives which are normal to the boundary, got (D[1](H))(x, 0)

 (1)

 (2)

 (3)

 (4)

## How can I send data from excel to maple in excel?...

Hello,

I'm trying to send data from Excel to Maple through Excel Add-in.

However, in Excel,  I can only create variables. I want to send these variables to Maple and access them in Maple.

I want to know how to send variables from Excel to Maple.

Thank you!

## Shift+Enter invalid in 2d input...

Hi, I have a procedure written on several lines on a worskheet in 2D. When I type shift+enter to insert a blank line at some place of the procedure, the blank line does not appear at the cursor, but at the end of the procedure.
This is annoying since I cannot add a line to make "readable" modifications in the procedure !

I am on LInux Fedora 26. Many thanks.

## Properties of functions and simplifications...

Hi all,

I'am new to Maple and have some very basic questions. Maple can figure out that the exponential is positive.

> is(exp(x)>0) assuming x::real;
true

But it cannot figure out that a sum of exponentials is positive too:

> is(sum(exp(x*i),i=1..N)>0) assuming x::real,N::integer,N>1;
FAIL

What am I doing wrong? Moreover, what's wrong with this statement:

> is(ln(a+b)-ln(a)=ln(1-b/a)) assuming a::real,b::real,a>0,b>0,b<a;
FAIL

Finally, is there a way to declare a generic fuction f and assume that its image is always a positive real?

Thanks

Nico

## How to prove a quadratic expression is positive...

Hi everybody,

Here is a very simple question :
I have two strictly positive reals r and s and the expression
a := r^2 – r*s + s^2
Obviously a is strictly positive too

But I’m not capable to prove this with Maple (2016, Windows 7), even after having declared r and s both positive which, I thought, would have prevent Maple to consider them as potentially complex.
I failed in using "is", "coulditbe", "use RealDomain ..." and so on

Totally distressed I tried to do something simpler ... which ultimately appeared still more upsetting

restart:
assume(r, ‘positive’):
coulditbe(r, ‘complex’) ;    # the answer is (surprising ?) true

I'm probably doing some mistake bigger than I but I don't find it !

Could you be compassionate enough to enlight me !!!

## finding standard error of the cofficents in nonlin...

hello i was wondering is it possible to find the standard error of the coefficents in the NonlinearFit function within the statistical package? or in any other package for that matter?

does someone have angorithme that can find it?

i've searched around the internet and it does not seem like maple employs this function at all?

hope someone can help:)

## How do I speed up the Evaluation of Elliptical fun...

Question:- can the procedure given below called "epi" be speeded up by compiling/ using evhf.If so how? My paple code is at the bottom.

First some background information.

Recently I ran into a difference in usage of a couple of elliptical functions between Maple and Mathematica.  This first case concerned EllipticalPi. The author of the blog kindly wrote  a Maple procedure to produce the same results as Mathematica’s  usage of ElllipticalPi.

I tested the basic integral that produces the EllipticPi    Ell := int(1/(1-nu*JacobiSN(t, k)^2), t)  answer      Ell := EllipticPi(JacobiSN(t, k), nu, k). They do not produce the same outcome. Plots are in the document .  They agree in one quarter only.

I then ran into a difference in usage of  EllipticF. This time I was able to get to same outcome myself using Maple’s help.

“It is worth noting the difference between the Legendre normal form of the Incomplete Elliptic integral of the first kind (see A&S 17.2.7), in Maple represented by EllipticF(z,k) but for the splitting of the square root in the denominator of the integrand (see definition lines above), and the normal trigonometric form of this elliptic integral (see A&S 17.2.6), in Maple represented by the InverseJacobiAM function
InverseJacobiAM(phi,k);

That worked fine.

There is no mention in the help for usage implementation of EllipticPi as opposed to different usages as there is with EllipticF. I do not know if there is a way in Maple of achieving the same enactment as Mathematica in this case, without the Procedure I  was  given.

Elliptic Pi in Mathematica and Maple

We use EllipticPi when we write exact solutions of rotation of a free asymmetric top. While solving Euler’s equations for angular velocity or angular momentum in the body frame we need Jacobi elliptic functions solving the differential equation for the attitude matrix involves EllipticPi function. As I have explained it in Taming the T-handle continued we need the integral

(1)

In Mathematica this is easily implemented as

(2)

However, as pointed out by Rowan in a comment to Taming the T-handle continued , the same formula does not work with Maple.

While the documentations of both Mathematica and Maple contain links to Abramowitz and Stegun Handbook of Mathematical Functions, they use different definitions. Here is what concerns us, from p. 590 of the 10th printing:

What we need is 17.2.16, while Maple is using 17.2.14. To convert we need to set but such a conversion is possible only in the domain where can be inverted. We can do it easily for sufficiently small values of but not necessarily for values that contain several quarter-periods.

The following Maple procedure does the job:

epi := proc (t::float, nu::float, k::float) local t2, n, dt, ep0, res; ep0 := EllipticPi(nu, k); t2 := EllipticK(k); n := floor(t/t2); dt := t-t2*n; if type(n, even) then res := Re(n*ep0+EllipticPi(JacobiSN(dt, k), nu, k)) else res := Re((1+n)*ep0-EllipticPi(JacobiSN(t2-dt, k), nu, k)) end if end proc

HAs an example here is the Maple plot for nu=-3, k=0.9:
plot(('epi')(t, -3.0, .9), t = -20 .. 20)

And here is the corresponding Mathematica plot:

The function epi(t,nu, k) defined above for Maple gives now the same result as EllipticPi(nu,JacobiAM(t,k^2),k^2) in Mathematica.

restart;
epi := proc (t, nu, k) local t2, n, dt, ep0, res; ep0 := EllipticPi(nu, k); t2 := EllipticK(k); n := floor(t/t2); dt := t-t2*n; if type(n, even) then res := Re(n*ep0+EllipticPi(JacobiSN(dt, k), nu, k)) else res := Re((1+n)*ep0-EllipticPi(JacobiSN(t2-dt, k), nu, k)) end if end proc;

Ell := int(1/(1-nu*JacobiSN(t, k)^2), t);
Ell := EllipticPi(JacobiSN(t, k), nu, k)
k := .9;
k := 0.9
nu := -3;
nu := -3
plot([epi(t, nu, k), Ell], t = -8 .. 20);

## fitting data on ode system...

hello agian i have the following problem. i need to fit data using the following model:

ode_sub := diff(S(t), t) = -k1*S(t)-S(t)/T1_s;

ode_P1 := diff(P1(t), t) = 2*k1*S(t)-k2*(P1(t)-P2(t)/keq)-P1(t)/T1_p1;

ode_P2 := diff(P2(t), t) = -k2*(-keq*P1(t)+P2(t))/keq-k4*P2(t)-P2(t)/T1_p2;

ode_P2e := diff(P2_e(t), t) = k4*P2(t)-P2_e(t)/T1_p2_e;

ode_system := ode_sub, ode_P1, ode_P2, ode_P2e

known paramters: s0 := 10000; k2 := 1000; T1_s := 14; T1_p2_e := 35; T1_p2 := T1_p1

initial conditions: init := S(0) = s0, P1(0) = 0, P2(0) = 0, P2_e(0) = 0

using the following data the fitting is fine:

T := [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 64, 66, 68, 70, 72, 74, 76, 78, 80, 82, 84, 86, 88, 90, 92, 94, 96, 98, 100]

data_s := [10000.00001377746, 7880.718836217512, 6210.572917625112, 4894.377814704496, 3857.121616806618, 3039.689036293312, 2395.493468824973, 1887.821030424934, 1487.738721378254, 1172.445015238064, 923.970970533911, 728.1555148262845, 573.8389058920359, 452.2263410725434, 356.3868133709972, 280.8584641247961, 221.3366863178145, 174.4291648589732, 137.4627967029932, 108.3305246342751, 85.37230370803576, 67.2794974831867, 53.0210755540487, 41.78436556408739, 32.92915589156605, 25.95049906181449, 20.45092590987601, 16.11678944747619, 12.70115104646253, 10.009439492356, 7.888058666357939, 6.216464754698855, 4.899036441885205, 3.860793948052506, 3.042584031379331, 2.397778731385364, 1.889601342133927, 1.489145259940784, 1.173591417634974, .9248694992255094, .7288443588090404, .5743879613705721, .4526024635132348, .3567687570029392, .2810508404011898, .2215650452813882, .174603181767829, .1376243372528356, .1084232889842753, 0.8542884707952822e-1, 0.6738282660463157e-1]

data_p1 := [0.1194335334401124e-4, 244.8746046039949, 374.8721199398692, 430.5392805383767, 439.6598364813143, 421.0353424914179, 387.1842556288343, 346.2646897222593, 303.4377508746471, 261.8283447091155, 223.1996547160051, 188.4213144493491, 157.7924449350029, 131.2622073983344, 108.5771635112278, 89.37951009190863, 73.26979150957087, 59.84578653950572, 48.72563358658898, 39.56010490461378, 32.03855466968536, 25.88922670933322, 20.87834763772145, 16.80708274458702, 13.50774122768974, 10.84014148654258, 8.687656394505874, 6.954093898245485, 5.560224107929433, 4.441209458524726, 3.544128529596104, 2.825755811619965, 2.251247757181308, 1.792233305651086, 1.425861347838012, 1.133566009019768, .90081361320016, .7153496336919163, .5678861241754847, .4505952916932289, .3572989037753538, .2832489239941939, .2244289868248577, .1778450590752305, .1408633578784151, .1114667192753896, 0.8826814044702111e-1, 0.6979315954603076e-1, 0.5526502783606788e-1, 0.4370354298880999e-1, 0.3456334307662573e-1]

data_p2_p2e := [-0.1821397630630296e-4, 1000.40572909871, 1568.064904416198, 1848.900129881268, 1944.147939710225, 1923.352973299286, 1833.705314342611, 1706.726235937363, 1563.036042115902, 1415.741121363331, 1272.825952816517, 1138.833091575137, 1016.03557293539, 905.2470623752856, 806.3754051707843, 718.7979384094563, 641.6091822001032, 573.7825966275556, 514.2718125966452, 462.0710416682647, 416.2491499591012, 375.9656328260581, 340.4748518743264, 309.124348579787, 281.3484612079911, 256.6599441255977, 234.6416876403397, 214.9378461256197, 197.2453333305823, 181.3059798685686, 166.90059218783, 153.8422174795751, 141.9717783871606, 131.1529759058513, 121.2692959115991, 112.2199678247825, 103.9187616370328, 96.29007054510998, 89.26877137303566, 82.79743967408479, 76.82586273439793, 71.30944943617081, 66.2085909515396, 61.48838150744805, 57.11714763242225, 53.06666224006544, 49.31130738219119, 45.82807853990728, 42.59597194910467, 39.59575450632008, 36.81013335261527]

the fitting is done the following way:

P1fu, P2fu, P2e_fu, Sfu := op(subs(res, [P1(t), P2(t), P2_e(t), S(t)]))

making residuals:

Q := proc (T1_p1, k1, keq, k4) local P1v, P2v, P2e_v, Sv, resid; option remember;

res(parameters = [T1_p1, k1, keq, k4]);

try P1v := ~[P1fu](T); P2v := ~[P2fu](T); P2e_v := ~[P2e_fu](T); Sv := ~[Sfu](T); resid := [P1v-data_p1, P2v+P2e_v-data_p2_p2e, Sv-data_s]; return [seq(seq(resid[i][j], i = 1 .. 3), j = 1 .. nops(T))] catch: return [1000000$3*nops(T)] end try end proc; q := [seq(subs(_nn = n, proc (T1_p1, k1, keq, k4) Q(args)[_nn] end proc), n = 1 .. 3*nops(T))]; finding inital point for the LSsolve: L := fsolve(q[2 .. 5], [10, 0.2e-1, 4, 4]) fitting the data with the intial point: solLS := Optimization:-LSSolve(q, initialpoint = L) this is all good however, when i used the following data it did not turn out so well (using the same approch as above): data_s := [96304.74567, 77385.03700, 62621.83067, 51239.94333, 42663.82367, 35084.74100, 28480.28367, 23066.01467, 18774.73700, 15179.13700, 12278.50767, 9937.652000, 8046.848333, 6521.242000, 5287.811667, 4277.779000, 3466.518333, 2835.467000, 2297.796333, 1861.249667, 1529.654000, 1235.353000, 999.6626667, 826.2343333, 667.9480000, 559.9230000, 449.2790000, 376.4860000, 289.1203333, 236.1483333] data_p1 := [0.86e-1, 3.904, 26.975, 31.719, 41.067, 46.779, 52.115, 43.101, 44.344, 41.094, 36.523, 27.742, 26.543, 28.062, 22.178, 21.303, 14.951, 17.871, 11.422, 12.051, 9.232, 6.817, 6.1, .717, 1.215, 6.146, .772, .375, 2.595, .518] data_p2_p2e := [-3.024, 22.238, 61.731, 103.816, 132.695, 159.069, 167.302, 160.188, 158.398, 152.943, 146.745, 135.22, 132.145, 120.413, 107.864, 95.339, 90.775, 81.828, 71.065, 70.475, 62.872, 49.955, 40.858, 42.938, 41.311, 35.583, 31.573, 29.841, 29.558, 21.762] the known parameters is the case are: s0 := 96304.74567; k2 := 10^5; T1_s := 14; T1_p2_e := 35; T1_p2 := T1_p1 additionally the following fuction affects the solution: K:=t->cos((1/180)*beta*Pi)^(t/Tr) i included this by doing this: P1fu_K := proc (t) options operator, arrow; P1fu(t)*K(t) end proc; P2fu_K := proc (t) options operator, arrow; P2fu(t)*K(t) end proc; P2e_fu_K := proc (t) options operator, arrow; P2e_fu(t)*K(t) end proc; Sfu_K := proc (t) options operator, arrow; Sfu(t)*K(t) end proc resulting in the following residuals: Q := proc (T1_p1, k1, keq, k4) local P1v, P2v, P2e_v, Sv, resid; option remember; res(parameters = [T1_p1, k1, keq, k4]); try P1v := ~[P1fu_K](T); P2v := ~[P2fu_K](T); P2e_v := ~[P2e_fu_K](T); Sv := ~[Sfu_K](T); resid := [P1v-data_p1, P2v+P2e_v-data_p2_p2e, Sv-data_s]; return [seq(seq(resid[i][j], i = 1 .. 3), j = 1 .. nops(T))] catch: return [1000000$3*nops(T)] end try end proc;
q := [seq(subs(_nn = n, proc (T1_p1, k1, keq, k4) Q(args)[_nn] end proc), n = 1 .. 3*nops(T))];

i think my problem is that the inital point in this case is not known. all i know is that all the fitted parameters should be positive and that k1<1, k4>10, keq>1 and T1_p1>100 (more i do not know) - is there a way to determin the inital point without guessing?

i also know the results, which should be close to these values:

k1=0.000438, k4=0.0385, keq=2.7385 and T1_p1=36.8 the output fit should look something like this

where the red curve is (PP2(t)+PP2_e(t))*K(t)

the blue cure is PP1(t)*K(t)

anyone able to help - i've tried for 2 days now. it might be that  ode_P1 := diff(P1(t), t) = 2*k1*S(t)-k2*(P1(t)-P2(t)/keq)-P1(t)/T1_p1 should be changed into ode_P1 := diff(P1(t), t) = k1*S(t)-k2*(P1(t)-P2(t)/keq)-P1(t)/T1_p1;

i've tried this but it didnt seem to do much

anyone able to help?:)

NB stiff=true can be used within the dsolve to speed up the process if needed:)