Advanced Search

MaplePrimes Questions

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

Asked by: illuminates 65
physics
April 05 2017
1  1

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])

Maple gives answer:

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))

Numerical integration...

Asked by: tazatel 50
eval add
April 05 2017
0  2

I have a derivative distribution function f that is defined:

f:=piecewise(x < .114e-1,0.,x < .129e-1,0.,x < .147e-1,0.,x < .167e-1,0.,x < .189e-1,0.,x < .215e-1,0.,x < .244e-1,0.,x < .278e-1,0.,x < .315e-1,0.,x < .358e-1,0.,x < .407e-1,0.,x < .463e-1,0.,x < .526e-1,0.,x < .597e-1,0.,x < .679e-1,0.,x < .771e-1,0.,x < .876e-1,0.,x < .995e-1,0.,x < .113,0.,x < .128,0.,x < .146,0.,x < .166,0.,x < .188,0.,x < .214,0.,x < .243,0.,x < .276,0.,x < .314,0.,x < .357,0.,x < .405,0.,x < .46,0.,x < .523,0.,x < .594,0.,x < .675,-.8800000000+1.481481481*x,x < .767,-.3935869565+.7608695652*x,x < .872,-.3213333333+.6666666667*x,x < .991,-.2529411765+.5882352941*x,x < 1.13,-.1690647482+.5035971223*x,x < 1.28,-.2026666667+.5333333333*x,x < 1.45,-.1223529412+.4705882353*x,x < 1.65,-.1650000000+.5000000000*x,x < 1.88,-.2726086957+.5652173913*x,x < 2.13,-.5636000000+.7200000000*x,x < 2.42,-.8662068966+.8620689655*x,x < 2.75,-1.200000000+1.000000000*x,x < 3.12,-1.645945946+1.162162162*x,x < 3.55,-1.865581395+1.232558140*x,x < 4.03,-2.075416667+1.291666667*x,x < 4.58,-1.925818182+1.254545455*x,x < 5.21,-1.559682540+1.174603175*x,x < 5.92,-.7967605634+1.028169014*x,x < 6.72,.3320000000+.8375000000*x,x < 7.64,1.942608696+.5978260870*x,x < 8.68,3.791923077+.3557692308*x,x < 9.86,5.850169492+.1186440678*x,x < 11.2,7.902985075-.8955223881e-1*x,x < 12.7,9.737333333-.2533333333*x,x < 14.5,10.82388889-.3388888889*x,x < 16.4,11.78631579-.4052631579*x,x < 18.7,11.34347826-.3782608696*x,x < 21.2,10.85240000-.3520000000*x,x < 24.1,9.311379310-.2793103448*x,x < 27.4,7.619090909-.2090909091*x,x < 31.1,5.814864865-.1432432432*x,x < 35.3,4.025714286-.8571428571e-1*x,x < 40.1,2.544375000-.4375000000e-1*x,x < 45.6,1.519090909-.1818181818e-1*x,x < 51.8,.8370967742-.3225806452e-2*x,x < 58.9,.6700000000,x < 66.9,.5963750000+.1250000000e-2*x,x < 76,.9005494505-.3296703297e-2*x,x < 86.4,1.161538462-.6730769231e-2*x,x < 98.1,1.318461538-.8547008547e-2*x,x < 111,1.544651163-.1085271318e-1*x,x < 127,1.241875000-.8125000000e-2*x,x < 144,1.106470588-.7058823529e-2*x,x < 163,.7721052632-.4736842105e-2*x,x < 186,0,x < 211,0,x < 240,0,x < 272,0,x < 310,0,x < 352,0,x < 400,0,x < 454,0,x < 516,0,x < 586,0,x < 666,0,x < 756,0,x < 859,0,x < 976,0,x < 1110,0,x < 1260,0,x < 1430,0,x < 1630,0,x < 1850,0,x < 2100,0,x < 2390,0,x < 2710,0,0):

I want to do a simple numerical integration of that function - that means I have k points of independent variable x and together with dependent variable points it looks like:

I was trying to resolve this in Maple:

f_list:=convert(f,list):

integralbodyx:=[seq(0+((200-0)/100)*i,i=1..100)]:

integralbodyy:=[seq(sum((eval(f_list[0+((200-0)/100)*j],x=(0+((200-0)/100)*j)))*2,j=1..i),i=1..100)]:

but I get an error 

Error, (in limit/mrv/limsimpl) too many levels of recursion

Can you help me with the syntax? Thank you.
 

I need help for these examples...

Asked by: kmasoud7787 10
homework
April 05 2017
1  0

Need numerical solution of PDE...

Asked by: Muhammad Usman 245
numeric pdsolve
April 05 2017
0  3

Dears 

Hope you would be fine. I want to solve the following PDEs by numerically for v[nf]=alpha[nf]=Ec=mu[nf]=C=1 and Pr=6.2

Eq1 := diff(u(x, t), t) = v[nf]*(diff(u(x, t), x, x));

Eq2 := diff(u(x, t), t) = alpha[nf]*(diff(theta(x, t), x, x))/Pr+Ec*mu[nf]*C*(diff(u(x, t), x))^2;

ICs := u(x, 0) = 0, theta(x, 0);

BCs := u(0, t) = 1, theta(0, t) = 1, u(10, t) = 0, theta(10, t) = 0;

and find the values of (diff(u(0, t), x))/(1-phi)^2.5 for different values of phi. Thanks in advace 

With my best regards and sincerely.

Muhammad Usman

School of Mathematical Sciences 
Peking University, Beijing, China

How can I wite a Maple code that can handle severa...

Asked by: owoicholohi 5
do-loop
April 04 2017
1  2

Please can someone help me with looping involving 4 variables such as S1 S2 S3 S4 from a series

Sn+1 =f(f-a)+u+Vn+w

Vn+1 =wc+f a-An+ u  if An = Un+Vn

Emacs Maple Mode...

Asked by: Thomas Dean 327
April 04 2017
0  2

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...

Asked by: Thomas Dean 327
plot
April 04 2017
0  1

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

plot(asdf);

actually plots y=x!

How to create a hyperplane which perpendicular to ...

Asked by: LeeHoYeung 535
April 04 2017
0  0

How to create a hyperplane which perpendicular to groebner basis

How do I use IsIsomorphic on large amount of small...

Asked by: tomasmedici 15
graph matrix isomorphism
April 04 2017
3  31

Hello,

I am trying to write a procedure to see, which (di)-graphs are isomorphic (here represented by there 3*3 adjecency-matrices). When I try the procedure for all 3*3-matrices with entries in {0,1} (there are 512 of them), I get the following error:
"Error, (in GraphTheory:-IsIsomorphic) invalid subscript selector"

Can you possibly say, what I am doing wrong? My code is the following:

getIso3 := proc(liste)
  local i,k,M1,c,d:
  c := 0:
  M1 := [[liste[1]]]:
  for i from 2 to numelems(liste) do
    for k from 1 to numelems(M1) do
        if IsIsomorphic(Digraph([a,b,c], liste[i]),Digraph([a,b,c],M1[k][1])) then
          M1[k] := [op(M1[k]), liste[i]]:
          c := 1:
        end if:
    end do:
    if c=0 then
      M1 := [op(M1), [liste[i]]]:
    else
      c := 0:
    end if:   end do: 
  return(M1):
end proc:

 

My input is a list (JJ) of 512 3*3 matrices constructed the following way :

all9Perm := proc(list)
  local P,i,m,n,A:
  P := list:
  for i from 0 to 9 do
    m:= i:  n:= 9-i:
    A := combinat:-permute([1$n, 0$m]):
    P := [op(P), op(1..numelems(A),A)]:
  od:
  return(P):
end proc:
K := []:
L := all9Perm(K):
listoflistsToListofmatrices := proc(liste)
  local M,i:
  M := []:
  for i from 1 to numelems(liste) do
      M := [op(M), Matrix([
          [ liste[i][1] , liste[i][2] , liste[i][3] ],
          [ liste[i][4] , liste[i][5] , liste[i][6] ],
          [ liste[i][7] , liste[i][8] , liste[i][9] ]])
          ]:
  end do:
  return(M):
end proc:
JJ := listoflistsToListofmatrices(L):

 

When I run this procedure on some of the 512 matrices it does work, but it crashes somewhere around matrix 350. I have try so split the list of the 512 matrices, and I am able to run the procedure on these splits, but this is very inconvenient :-)

I hope you can help me. Also if this can be done in an easier way - I am new to programming and recieve help with a smile.

Yours, Tomas.

Error, (in Groebner:-NormalSet) The case of non-ze...

Asked by: LeeHoYeung 535
groebner
April 04 2017
0  3
tord := plex(x, y, z);
G := Basis([hello1, hello2, hello3], tord);
ns, rv := NormalSet(G, tord);
Error, (in Groebner:-NormalSet) The case of non-zero-dimensional varieties is not handled
is this error due to version of maple?
which version do not have this error?
 

How to solve symbolically for like following coupl...

Asked by: student_1 305
April 04 2017
0  3

I wonder how to solve symbolically for like following coupled ODEs in Maple? 

On the other hand, I want to write the code for step by step solution of this problem. But I didn't find algorithm of the solution on the some books. Do you know some books including solving coupled ODEs ?

 

How to solve delay differential equation by method...

Asked by: rkdevkate 5
dsolve dde
April 04 2017
1  2

How to solve delay differential equation by method of steps in MAPLE software. 

How do I produce a sequence defined by a linear re...

Asked by: st104290 35
recurrence sequence
April 04 2017
2  5

In Mathematica, there is a command "LinearRecurrence" that lists a sequence generated by a linear recurrence.
For example, 
LinearRecurrence[{-3, 1}, {7, 2}, 10]
will produce the first ten numbers in the sequence defined by a0=7, a1=2; a_n = -3 an-1 + an-2.

Is there anything similar in Maple? 
If not, what would be the easiest way to do the same?

UPDATE: Two answers show how to do this example with rsolve, but it seems that this solution does not help in general for recurrence of high order.

Thank you.

How can I generate series of variables using Maple...

Asked by: Elisha 50
do-loop
April 04 2017
0  0


 

> 

restart
> 

with(student)
> 

``
> 

G := S(t)*L(t)

S(t)*L(t)

 (1)
> 

m := 10
> 

S[lambda] := sum(S[b]*lambda^b, b = 0 .. m); L[lambda] := sum(L[b]*lambda^b, b = 0 .. m); G[lambda] := subs(S(t) = S[lambda], G); G[lambda] := subs(L(t) = L[lambda], G[lambda]); G := G[lambda]; s := expand(G, lambda); ft := unapply(s, lambda); for i from 0 while i <= m do A1[i] := ((D@@i)(ft))(0)/factorial(i); print(A[i] = A1[i]) end do

A[0] = S[0]*L[0]

  

A[1] = L[0]*S[1]+L[1]*S[0]

  

A[2] = L[0]*S[2]+L[1]*S[1]+L[2]*S[0]

  

A[3] = L[0]*S[3]+L[1]*S[2]+L[2]*S[1]+L[3]*S[0]

  

A[4] = L[0]*S[4]+L[1]*S[3]+L[2]*S[2]+L[3]*S[1]+L[4]*S[0]

  

A[5] = L[0]*S[5]+L[1]*S[4]+L[2]*S[3]+L[3]*S[2]+L[4]*S[1]+L[5]*S[0]

  

A[6] = L[0]*S[6]+L[1]*S[5]+L[2]*S[4]+L[3]*S[3]+L[4]*S[2]+L[5]*S[1]+L[6]*S[0]

  

A[7] = S[5]*L[2]+S[2]*L[5]+S[1]*L[6]+S[0]*L[7]+S[4]*L[3]+S[3]*L[4]+S[6]*L[1]+S[7]*L[0]

  

A[8] = S[2]*L[6]+S[0]*L[8]+S[8]*L[0]+S[4]*L[4]+S[1]*L[7]+S[3]*L[5]+S[6]*L[2]+S[5]*L[3]+S[7]*L[1]

  

A[9] = S[6]*L[3]+S[9]*L[0]+S[4]*L[5]+S[2]*L[7]+S[3]*L[6]+S[8]*L[1]+S[0]*L[9]+S[1]*L[8]+S[7]*L[2]+S[5]*L[4]

  

A[10] = S[0]*L[10]+S[10]*L[0]+S[1]*L[9]+S[6]*L[4]+S[7]*L[3]+S[3]*L[7]+S[2]*L[8]+S[4]*L[6]+S[9]*L[1]+S[5]*L[5]+S[8]*L[2]

 (2)

s[n+1] := (1-f)*alpha*(int(v__n, t = 0 .. t))-beta*c*(int(A__n, t = 0 .. t))-(`&theta;__1`+a+Pi)*(int(s__n, t = 0 .. t))

v[n+1] := `&theta;__1`*(int(s__n, t = 0 .. t))-((1-f)*alpha+f*`&theta;__2`+a+Pi)*(int(v__n, t = 0 .. t))

e[n+1] := `&beta;c`*(int(A__n, t = 0 .. t))-(delta+a+Pi)*(int(e__n, t = 0 .. t))

i[n+1] := delta*(int(e__n, t = 0 .. t))-(eta+a+Pi)*(int(i__n, t = 0 .. t))

r[n+1] := eta*(int(i__n, t = 0 .. t))+`f&theta;__2`*(int(v__n, t = 0 .. t))-(a+Pi)*(int(r__n, t = 0 .. t))

for n from 0 to 4 do end do

``

(1-f)*alpha*v__n*t-beta*c*s__0*i__0*t-(`&theta;__1`+a+Pi)*s__n*t

 (3)

"for n:=0, n=1, n=2, n=3, n=4: `s__1, ``s__2__`, `s__3, ``s__4`, `i__1`, `i__2, ``i__3`, `i__4`, `r__1`, `r__2`, `r__3, ``r__4`:"

Error, unterminated for loop

"for n:=0, n=1, n=2, n=3, n=4: `s__1, ` `s__2__`, `s__3, __3__, s__4`, `i__1__`, `i__2, ` `i__3`, `i__4`, `r__1`, `r__2`, `r__3, ` `r__4`:"

  

"and s(t)=`s__1_`+`s__2 __`+ `s__3, `+`s__4`,"

i(t) = i__1+`i__2, `+i__3+i__4, r(t) = r__1+r__2+r__3_+r__4

"but A[n]:=(1)/(n!) [((&DifferentialD;)^(n))/(&DifferentialD; lambda^(n)) ((&sum;)`s__k`lambda^(k))((&sum;)`i__k`lambda^(k)) ]() ? ()|() ? (lambda=0)"

Error, missing numerator

Typesetting:-mambiguous(Typesetting:-mrow(Typesetting:-mi("but", fontstyle = "italic", mathvariant = "italic"), Typesetting:-mo("&InvisibleTimes;", accent = "false", fence = "false", largeop = "false", lspace = "0.0em", mathvariant = "normal", movablelimits = "false", rspace = "0.0em", separator = "false", stretchy = "false", symmetric = "false"), Typesetting:-msub(Typesetting:-mi("A", fontstyle = "italic", mathvariant = "italic"), Typesetting:-mi("n", fontstyle = "italic", mathcolor = "#c800c8", mathvariant = "italic", placeholder = "true"), subscriptshift = "0"), Typesetting:-mo("&Assign;", accent = "false", fence = "false", largeop = "false", lspace = "0.2777778em", mathvariant = "normal", movablelimits = "false", rspace = "0.2777778em", separator = "false", stretchy = "false", symmetric = "false"), Typesetting:-mfrac(Typesetting:-mn("1", mathvariant = "normal"), Typesetting:-mrow(Typesetting:-mi("n", fontstyle = "italic", mathvariant = "italic"), Typesetting:-mo("&excl;", accent = "false", fence = "false", largeop = "false", lspace = "0.1111111em", mathvariant = "normal", movablelimits = "false", rspace = "0.1111111em", separator = "false", stretchy = "false", symmetric = "false")), bevelled = "false", denomalign = "center", linethickness = "1", numalign = "center"), Typesetting:-mo("&InvisibleTimes;", accent = "false", fence = "false", largeop = "false", lspace = "0.0em", mathvariant = "normal", movablelimits = "false", rspace = "0.0em", separator = "false", stretchy = "false", symmetric = "false"), Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mfrac(Typesetting:-msup(Typesetting:-mo("&DifferentialD;", accent = "unset", fence = "unset", largeop = "unset", lspace = "0.0em", mathvariant = "normal", movablelimits = "unset", rspace = "0.0em", separator = "unset", stretchy = "unset", symmetric = "unset"), Typesetting:-mi("n", fontstyle = "italic", mathvariant = "italic"), superscriptshift = "0"), Typesetting:-mrow(Typesetting:-mo("&DifferentialD;", accent = "unset", fence = "unset", largeop = "unset", lspace = "0.0em", mathvariant = "normal", movablelimits = "unset", rspace = "0.0em", separator = "unset", stretchy = "unset", symmetric = "unset"), Typesetting:-mo("&InvisibleTimes;", accent = "false", fence = "false", largeop = "false", lspace = "0.0em", mathvariant = "normal", movablelimits = "false", rspace = "0.0em", separator = "false", stretchy = "false", symmetric = "false"), Typesetting:-msup(Typesetting:-mi("&lambda;", fontstyle = "normal", mathvariant = "normal"), Typesetting:-mi("n", fontstyle = "italic", mathvariant = "italic"), superscriptshift = "0")), bevelled = "false", denomalign = "center", linethickness = "1", numalign = "center"), Typesetting:-mo("&InvisibleTimes;", accent = "false", fence = "false", largeop = "false", lspace = "0.0em", mathvariant = "normal", movablelimits = "false", rspace = "0.0em", separator = "false", stretchy = "false", symmetric = "false"), Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-munderover(Typesetting:-mo("&sum;", accent = "unset", fence = "unset", largeop = "true", lspace = "0.0em", mathvariant = "normal", movablelimits = "true", rspace = "0.1666667em", separator = "unset", stretchy = "true", symmetric = "unset"), Typesetting:-mrow(Typesetting:-mi("k", fontstyle = "italic", mathvariant = "italic"), Typesetting:-mo("&equals;", accent = "false", fence = "false", largeop = "false", lspace = "0.2777778em", mathvariant = "normal", movablelimits = "false", rspace = "0.2777778em", separator = "false", stretchy = "false", symmetric = "false"), Typesetting:-mn("0", mathvariant = "normal")), Typesetting:-mn("4", mathvariant = "normal"), accent = "false", accentunder = "false"), Typesetting:-mi("`s__k`"), Typesetting:-msup(Typesetting:-mi("&lambda;", fontstyle = "normal", mathvariant = "normal"), Typesetting:-mi("k", fontstyle = "italic", mathvariant = "italic"), superscriptshift = "0")), mathvariant = "normal"), Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-munderover(Typesetting:-mo("&sum;", accent = "unset", fence = "unset", largeop = "true", lspace = "0.0em", mathvariant = "normal", movablelimits = "true", rspace = "0.1666667em", separator = "unset", stretchy = "true", symmetric = "unset"), Typesetting:-mrow(Typesetting:-mi("k", fontstyle = "italic", mathvariant = "italic"), Typesetting:-mo("&equals;", accent = "false", fence = "false", largeop = "false", lspace = "0.2777778em", mathvariant = "normal", movablelimits = "false", rspace = "0.2777778em", separator = "false", stretchy = "false", symmetric = "false"), Typesetting:-mn("0", mathvariant = "normal")), Typesetting:-mn("4", mathvariant = "normal"), accent = "false", accentunder = "false"), Typesetting:-mi("`i__k`"), Typesetting:-msup(Typesetting:-mi("&lambda;", fontstyle = "normal", mathvariant = "normal"), Typesetting:-mi("k", fontstyle = "italic", mathvariant = "italic"), superscriptshift = "0")), mathvariant = "normal"), Typesetting:-mo("&InvisibleTimes;", accent = "false", fence = "false", largeop = "false", lspace = "0.0em", mathvariant = "normal", movablelimits = "false", rspace = "0.0em", separator = "false", stretchy = "false", symmetric = "false")), open = "&lsqb;", close = "&rsqb;", mathvariant = "normal"), Typesetting:-mfrac(Typesetting:-mambiguous(Typesetting:-merror("?"), Typesetting:-merror("missing numerator")), Typesetting:-mphantom(Typesetting:-mrow(Typesetting:-mi("x", fontstyle = "italic", mathvariant = "italic"), Typesetting:-mo("&equals;", accent = "unset", fence = "unset", largeop = "unset", lspace = "0.2777778em", mathvariant = "normal", movablelimits = "unset", rspace = "0.2777778em", separator = "unset", stretchy = "unset", symmetric = "unset"), Typesetting:-mi("a", fontstyle = "italic", mathvariant = "italic")), constraints = "height-only"), bevelled = "false", denomalign = "center", linethickness = "0", numalign = "center"), Typesetting:-mo("&verbar;", accent = "unset", fence = "unset", largeop = "unset", lspace = "0.1111111em", mathvariant = "normal", movablelimits = "unset", rspace = "0.1111111em", separator = "unset", stretchy = "true", symmetric = "unset"), Typesetting:-mfrac(Typesetting:-mphantom(Typesetting:-mi("f(x)", fontstyle = "italic", mathvariant = "italic"), constraints = "height-only"), Typesetting:-mrow(Typesetting:-mi("&lambda;", fontstyle = "normal", mathvariant = "normal"), Typesetting:-mo("&equals;", accent = "false", fence = "false", largeop = "false", lspace = "0.2777778em", mathvariant = "normal", movablelimits = "false", rspace = "0.2777778em", separator = "false", stretchy = "false", symmetric = "false"), Typesetting:-mn("0", mathvariant = "normal")), bevelled = "false", denomalign = "center", linethickness = "0", numalign = "center")))

  

NULLI have tried to do the above but got the error messages. I tried changing the L in that computation above to i to enable me get the desired result but i couldn't.

> Please how can I generate the values S1,S2,S3,S4; i1,i2,i3,i4 and r1,r2,r3,r4 and finally do

S(t)=S1 + S2 + S3 + S4; i(t)= i1+i2+i3+i4 ; and r(t)=r1,r2,r3,r4  using Maple?

> 

NULL
> 

> 

> 

``
> 

``


 

Download Adomian.Elisha2.mw

How to plot a prettier bifurcation diagram using M...

Asked by: ZWang 25
imagetools iterative-map
April 04 2017
0  5

Dear All,

The following code plots the bifurcation diagram for a three-dimensional continuous dynamical system as a variable Re varies. However, the resulting plot (by pointplot command) is rather ugly, comparing with other bifurcation diagrams, see attached. Could anyone point me out how to improve its looking?

``

> 

restart:
with(plots): with(DEtools): with(plottools):with(LinearAlgebra):
doSol:=proc( R )
             local t_start:=2500,
                   t_end:=5000,
                   b:=-3,
                   c:=3,
                   v1:=1,
                   f:=-4,
                   v2:=2.0,
                   omega:=0.1*sqrt(R),
                   epsilon:=evalf(1/R),
                   k:=0.1,
                   kH:=(c+1)/(b-1),
                   sys, evs, w, M, T, i, tt, solA, N_alpha,
                   r_center, z_center, Zshift, alpha, alpha1,
                   alpha2, del_alpha, m, Z, Rshift, RR;
             sys:=diff(u1(t),t)=v1*u1(t)-(omega+k*u2(t)^2)*u2(t)-(u1(t)^2+u2(t)^2-b*z(t)^2)*u1(t),
                  diff(u2(t),t)=(omega+k*u1(t)^2)*u1(t)+v1*u2(t)-(u1(t)^2+u2(t)^2-b*z(t)^2)*u2(t),
                  diff(z(t),t)=z(t)*(kH*v1+c*u1(t)^2+c*u2(t)^2+z(t)^2)+epsilon*z(t)*(v2+f*z(t)^4):
             solA:=dsolve({sys, u1(0)=0.6, u2(0)=0.6, z(0)=0.1},
                          {u1(t),u2(t),z(t)},
                          numeric, method=rkf45, maxfun=10^7,
                          events=[[[u1(t)=0, u2(t)>0], halt]]
                         );
             evs:=Array():
             evs(1,1..4):=Array([t,u1(t),u2(t),z(t)]);
             interface(warnlevel=0):

             for i from 2 do
                 w:=solA(t_end):
                 if   rhs(w[1])<t_end
                 then evs(i,1..4):=Array(map(rhs, w));
                      solA(eventclear);
                 else break;
                 fi
             od:

             interface(warnlevel=3):
             M:=DeleteRow(convert(evs,matrix),1):

             T:=M[..,1]:

             for i from 1 do
                 tt:=T[i]:
                 if   tt>=t_start
                 then break;
                 end if;
             end do:

             RR:=M[i..,3]: Z:=M[i..,4]:

             N_alpha:=numelems(RR):

             r_center:=add(RR[i],i=1..N_alpha)/N_alpha:
             z_center:=add(Z[i],i=1..N_alpha)/N_alpha:

             Zshift:=Z-~r_center: Rshift:=RR-~z_center:

             alpha:=(arctan~(Zshift/~Rshift)+(Pi/2)*sign~(Rshift))/~(2*Pi)+~0.5:

             return alpha;
  end proc:
> 

i_start:=125: i_end:=250: i_delta:=0.1:

M:=[seq( [j, doSol(j)], j=i_start..i_end, i_delta)]:
> 

S:=seq([Vector(numelems(M[i,2]),fill=M[i,1]),M[i,2]], i=1..numelems(M),1):
display(seq( pointplot(S[i],
             symbolsize=1,
             symbol=point) ,
             i=1..numelems(M)),
             axes=boxed,
             view=[i_start..i_end,0..1],
             size=[650,400],
             axesfont= ["TimesNewRoman", 16],
             labels=["Re",phi[n]],
             labelfont = ["TimesNewRoman", 16],
             labeldirections=[horizontal, vertical]);

``


 

Bifurcation_diagram.mw

Thank you.

Very kind wishes,

Wang Zhe

First 990 991 992 993 994 995 996 Last Page 992 of 2428