tomleslie

4487 Reputation

10 Badges

9 years, 98 days

MaplePrimes Activity


These are answers submitted by tomleslie

quite a lot of "spurious" charactersin the string to be converted - eg \[ ]. Why do these exist? Some sort of copy/paste issue maybe? No idea how the translator may interpret these characters.

Maybe also have a problem with capitalisation??

The attached seems to work as I would expect

restart:

with(MmaTranslator):
interface(rtablesize=10):

Gamma

(1)

MM:=Matrix(convert("{{1/2 (-2 kappa - 4 lambda), -2 Sqrt[2] g, 0}, {g/Sqrt[
  2], -Gamma - kappa/2 - lambda, -Sqrt[2] g}, {0,
  Sqrt[2] g, -2 Gamma}}", FromMma));
simplify~(MM);

Matrix(3, 3, {(1, 1) = -kappa-2*lambda, (1, 2) = -2*sqrt(2)*g, (1, 3) = 0, (2, 1) = g/sqrt(2), (2, 2) = -GAMMA-(1/2)*kappa-lambda, (2, 3) = -sqrt(2)*g, (3, 1) = 0, (3, 2) = sqrt(2)*g, (3, 3) = -2*GAMMA})

 

Matrix(%id = 18446744074370869238)

(2)

 

 

 


 

Download mmaToMaple.mw

In the attached, you will need to convert the filepath to something appropriate for your installation.

There are many options which can be added to the "pointplot()" command

restart;

interface(rtablesize=10):

M:= ExcelTools:-Import("C:/Users/TomLeslie/Desktop/poljska_mreza.xlsx"):
plots:-pointplot( M[..,1..2]);

 

 

Download xlPlot.mw

You have three independent variables (x, y, t). The help page for pdsolve/numeric states clearly (my emphasis)

PDEsys
-
single or set or list of time-dependent partial differential equations in two independent variables

So you will not be able to solve your system numerically (unless some kind of symmetry allows you to reduce it to two independent variables).

So numerical solution is (probably) a non-starter.

You *may* be able to solve the system analytically - but if you want anyone to work on this, then upload your worksheet using the big green up-arrow in the Mapleprimes toolbar. Life is too short fro anyone here to retype your worksheet from a picture!

There is no "valuesplit" option associated with "colorscheme" in Maple 18 (which you state you are using). In fact the "colorscheme" capability in Maple 18 seems to be very, very limited.

So where did you get the idea of using "valuesplit" with "colorscheme"???? Maybe you are using Maple 2018, not Maple 18???

The attached shows a few options for graph coloring with "colorscheme". I'm not really sure what coloring you are trying to achieve, so the contents of the attached are pretty much "guesses".

Interestingly, the color variation in the various graphs is not reproduced when they are rendereed on this site, so you will have to download  the worksheet to examine

  restart:

  with(VectorCalculus):
  kernelopts(version);
  interface(rtablesize=10):

`Maple 2019.0, X86 64 WINDOWS, Mar 9 2019, Build ID 1384062`

(1)

  a := 1.0:
  q := a*cos(b):
  p := a*sin(b):
  R1 := PositionVector([q, p], cartesian[x, y]):
  PlotPositionVector( R1,
                      b = 0 .. 2*Pi,
                      curveoptions=[ colorscheme=[ "valuesplit",
                                                   [ -1..-1/2="Red",
                                                     -1/2..1/2="Green",
                                                      1/2..1="Blue"
                                                   ]
                                                 ]
                                   ]
                    );
  PlotPositionVector( R1,
                      b = 0 .. 2*Pi,
                      curveoptions=[ colorscheme=[ "linear",
                                                   ["Red","Blue"]
                                                 ]
                                   ]
                    );
  PlotPositionVector( R1,
                      b = 0 .. 2*Pi,
                      curveoptions=[ colorscheme=[ "linear",
                                                   ["Red", "Yellow", "Blue"]
                                                 ]
                                   ]
                    );

 

 

 

 

 


 

Download colsch.mw

just pass the index as an extra parameter. After all does it really matter to ou whether you call

f[3](x, y)  or

f(x, y, 3) 

If you can accept the latter as a calling sequence, then the function definition becomes

f := (x, y, i) -> a[i]*x + b[i]*y;

As far as I can tell, pretty much everything is working - so what exactly is the problem. See the attached

  restart;

  with(plots):
  interface(rtablesize=10):
  c1:= 1: c2:= 1: nu:= 3: Ra:= 2: sc:= 1: K:= [black, red, green]:
  A:= c2/(nu*c1): R:= Ra: N:= [0, 1, 2]: pr:= 1: S:= sc: g:= 0.5: k:= 1:
  for j to nops(N) do
      sol1:= dsolve
             ( [ diff(f(eta), eta$3)
                 +
                 A*( 3/4*f(eta)*diff(f(eta), eta$2)
                     -
                     1/2*diff(f(eta), eta)^2
                   )
                 +
                 R*( theta(eta) + N[j]*phi(eta) )
                 -
                 k*( 2*diff(f(eta), eta)
                     +
                     c1*diff(f(eta), eta$2)*diff(f(eta), eta$2)
                     -
                     3*diff(f(eta), eta)*diff(f(eta), eta$4)
                   )
                 = 0,
                 diff(theta(eta), eta$2)/pr + 3/4*f(eta)*diff(theta(eta), eta) = 0,
                 diff(phi(eta), eta$2)/S + 3/4*f(eta)*diff(phi(eta), eta) = 0,
                 f(0) = 0, D(f)(0) = 1, D(f)(3) = 1, (D@@2)(f)(3) = 0, phi(0) = 1,
                 phi(3) = 0, theta(0) = D(theta)(0)/g + 1, theta(3) = 0
               ],
               numeric,
               method = bvp[middefer]
             );
      fplt[j]:= odeplot(sol1, [eta, f(eta)], color = K[j], axes = boxed);
      tplt[j]:= odeplot(sol1, [eta, theta(eta)], color = K[j], axes = boxed);
      fplt[j]:= odeplot(sol1, [eta, diff(f(eta), eta)], color = K[j], axes = boxed);
      phiplt[j]:= odeplot(sol1, [[eta, phi(eta)]], color = K[j], axes = boxed);
  end do:
  display([seq(fplt[j], j = 1 .. nops(N))]);
  display([seq(tplt[j], j = 1 .. nops(N))]);
  display([seq(phiplt[j], j = 1 .. nops(N))]);

 

 

 

#
# OP's second(?) problem
#
  sol2:= dsolve
         ( [ diff(f(eta), eta$3)+ f(eta)* diff(f(eta), eta$2)=0,
             diff(theta(eta), eta$2)+ f (eta)* diff(theta(eta), eta)=0,
             f(0) = 0, (D(f))(0) = 0, (D(f))(5) = 1,
             theta(0) = 1, theta(5) = 0
           ],
           numeric,
           method=bvp[midrich]
         ):
  odeplot( sol2,
           [ [eta, f(eta)],
             [eta, theta(eta)]
           ],
           eta=0..3,
           color=[red, green]
         );

 

 

Download odeProb.mw

using Maple2018.2 (Build ID 1362973) on Windows 7.

See the attached

restart;

interface(version);
with(TimeSeriesAnalysis):
esm2 := ExponentialSmoothingModel(seasonal = {A, M}, constraints = admissible);

`Standard Worksheet Interface, Maple 2018.2, Windows 7, November 16 2018 Build ID 1362973`

 

_m652943904

(1)

 

Download timeSer.mw

Not so much a singularity, more the fact that the system is inconsistent.

Each differential in 'odesys' is multiplied by the factor sqrt(t). Hence at t=0, all of these differential terms disappear, and you are left with a set of linear equations in S(0), M(0)., N(0)., U(0). Into these equations just substitute your boundary conditions, using

eval([eval(odesys, t=0)], [ICS])

which will return

[32.619 = 0, -5.19 = 0, -22.81 = 0, -0.30 = 0]

This is definitely not good!

see the attached

restart;

interface(rtablesize=10):

with(Groebner):
Jugdement := proc (F, P)
                   local f, N1, N2, i, j, k, q, LTF, LTPj, F__loc:=F;
                   N1 := nops(F__loc);
                   N2 := numelems(P);
                   for i to N1 do
                       LTF := LeadingTerm(F__loc, plex(w, t, s, x, y, z));
                       for j to N2 do
                           LTPj := LeadingTerm(P[j], plex(w, t, s, x, y, z));
                           f := divide(LTF[2], LTPj[2], 'q');
                           if   f = true
                           then k := j;
                                break
                           end if
                       end do;
                       if   f = true
                       then break
                       end if;
                       if   f = false
                       then F__loc := F__loc-LTF[1]*LTF[2]
                       end if
                   end do;
                   return k, f, q, P[k], LTF[1]
             end proc:
F := x^2+y^2+z^2;
P := {z, x*z, y*x};
Jugdement(F, P);

x^2+y^2+z^2

 

{z, x*z, y*x}

 

1, true, z, z, 1

(1)

 


 

Download formPar.mw

The GraphTheory() package does not permit Graphs with loops.

From the help page of ?GraphTheory (with my emphasis)

The GraphTheory package is a collection of routines for creating graphs, drawing graphs, manipulating graphs, and testing graphs for properties. The graphs are sets of vertices (nodes) connected by edges. The package supports both directed and undirected graphs but not multigraphs or graphs with loops (self-incident vertices). The edges in the graphs can be weighted or unweighted.

 

 

 - but you don't have to!

  1. The attached generates the desired plot from the original expresssion (ie using Cartesian coordinates)
  2. Then uses changecoords() to convert the expression to spherical polars,  and produces the "same" plot using the polar coordinates

  restart:

  interface(rtablesize=10):
  with(plots):

#
# Produce plot using cartesian coordinates
#
  expr1:=x^2+y^2/4+z^2/9=1:
  implicitplot3d( expr1,
                  x=-4..4,
                  y=-4..4,
                  z=-4..4,
                  axes=none,
                  style=surface,
                  scaling=constrained,
                  grid=[ 40, 40, 40]
               )

 

#
# Convert expression to spherical polars
#
  expr2:= :-changecoords( expr1,
                          [ x, y, z ],
                          spherical,
                          [ r, theta, phi ]
                        );
#
# Produce plot using spherical coordinates
#
  implicitplot3d( expr2,
                  r=0..4,
                  theta=0..2*Pi,
                  phi=0..Pi,
                  style=surface,
                  coords=spherical,
                  axes=none,
                  scaling=constrained,
                  style=surface,
                  grid=[40, 40, 40]
                );

r^2*sin(phi)^2*cos(theta)^2+(1/4)*r^2*sin(phi)^2*sin(theta)^2+(1/9)*r^2*cos(phi)^2 = 1

 

 

 


 

Download elliPlot.mw

The attached "just works" - and it isn't a version issue, becuase I've checked it in Maple 2018, Maple 2017 ,Maple 2016, Maple 2015, and Maple 18.

  restart;

  kernelopts(version);
  Physics:-Version();
  interface(rtablesize=10):

`Maple 2019.0, X86 64 WINDOWS, Mar 9 2019, Build ID 1384062`

 

"C:\Users\TomLeslie\maple\toolbox\2019\Physics Updates\lib\Physics Updates.maple", `2019, April 16, 8:7 hours, version in the MapleCloud: 344, version installed in this computer: 344.`

(1)

  PDE:=[ diff(f(x,xp),x)=-(1/2)*(L*xp+2*x)*kl,
         diff(f(x,xp),xp)=-(1/4)*kl*L*(L*xp+2*x)
       ];
  sol:=pdsolve(PDE);
  pdetest(sol[], PDE);

[diff(f(x, xp), x) = -(1/2)*(L*xp+2*x)*kl, diff(f(x, xp), xp) = -(1/4)*kl*L*(L*xp+2*x)]

 

{f(x, xp) = -(1/8)*(L*xp+2*x)^2*kl+_C1}

 

[0, 0]

(2)

 

 


 

Download pdprob.mw

which are fixed in the attached

#
# Restar should be in its own execution group
#
  restart;

#
# redefining built-in variables is really not
# advisable!
#
  local gamma;
  local I;
  local pi ;

gamma

 

I

 

Warning, The imaginary unit, I, has been renamed _I

 

I

 

pi

(1)

#
# Added this execution group so worksheet displays
# correctly on Mapleprimes website
#
  interface(rtablesize=10):

#
# Set up numerical values for all problem parameters
#
# Changed beta[*] to beta[star] in this execution group
# (and everywhere subsequently).
#
# rmoved excess comma character after final entry in params
#
  params:=[       psi=0.142,        mu[1]=0.112,      phi=0.4e-3,
                 mu[v]=0.002, beta[o]=0.081,  M[h]=10,
            omega=0.2e-2,     eta=0.5e-1, mu[e]=0.092,
                pi=0.598e-2,    beta[star]=.5,      eta=0.213
             
          ]:

#
# Define main function
#
  R:= sqrt((psi+mu[1]+phi)*(mu[1])^(2)*mu[v]*psi*beta[o]*(M[h])^(2)*(omega+mu[1]+eta)*mu[e]*pi*beta[star]/(psi+mu[1]+phi)*(omega+mu[1]+eta)*mu[v]*mu[e]*(mu[1])^2);

(mu[1]^4*mu[v]^2*psi*beta[o]*M[h]^2*(omega+mu[1]+eta)^2*mu[e]^2*pi*beta[star])^(1/2)

(2)

#
# Compute "all" derivatives and evaluate numerically.
#
# For the purposes of this calculation "all"
# derivatives, means the derivatives with respect to
# every variable returned by indets(R, name)
#
# Output a list of two element lists where each of
# the latter is
#
# [ varName,
#   eval( diff( R, varName), params )
# ]
#
 [ seq( [j, eval( diff( R, j), params )],j in indets(R, name))];

[[eta, 0.1353555737e-6], [omega, 0.1353555737e-6], [pi, 0.1856046329e-5], [psi, 0.7816307780e-7], [M[h], 0.2219831409e-8], [beta[o], 0.1370266302e-6], [beta[star], 0.2219831409e-7], [mu[1], 0.5317540395e-6], [mu[e], 0.2412860227e-6], [mu[v], 0.1109915705e-4]]

(3)

#
# Compute all "sensitivities" (where the sensitivity
# is as defined in Rouben Rostamian response to the
# OP's earlier post) and evaluate numerically.
#
# For the purposes of this calculation "all" sensitivities
# means the sensitivity with respect to every variable
# returned by indets(R, name)
#
# Output a list of two element lists where each of
# the latter is
#
# [ varName,
#   eval( varName*diff( R, varName)/R, params )
# ]
#
  seq( [j, eval( j*diff( R, j)/R, params )],j in indets(R, name));

[eta, .3048780488], [omega, 0.1219512195e-1], [pi, 1/2], [psi, 1/2], [M[h], 1], [beta[o], 1/2], [beta[star], 1/2], [mu[1], 2.682926826], [mu[e], 1], [mu[v], 1]

(4)

 

NULL


 

Download sens.mw

  1. You have used x both as the independent variable in all of the ODEs, but also as an index in a for-loop. The latter use will result on 5.2 being assigned to 'x' after the loop has executed. This will effectively substitute 5.2 everywhere you use 'x' in your ODEs. You really do not want this, so I have changed the offending loop index at the point highlighted in the attached
  2. Second problem is more subtle. Location is also highlighted in the attached. The loop which solves the ODEs for coefficients of p0, p1,p2,...etc. For p0 one gets an explicit expression fro c[0](x), whihc is good. However for p1, the relevant ODE contains c[0](x) (which will be correctly substituted from the previous loop iteration), c[1](x) and f(x). It is therefore not possible to get an explicit expression for c[1](x). Maple will provide a 'formal' solution for c[1](x) in terms of double integrals of the unknown f(x) - but I'm guessing that this is not what you want! It gets worse! This formal soultion for c[1](x) is used in the ODE for p2 which also contains the two unknowns c[2](x) and f(x). Maple will provide another 'formal' solution, but this contains the previous 'formal' solution for c[1](x) - so ones gets an even more complicated formal solution containing a formal solution. This situation gets worse at each loop iteration, and the expression complexity "explodes". I'm guessing thsat this is not what you want, but it is an inevitable consequence of the way the problem is formulated. Can only suggest that you check the logic of your algorithm very carefully

NULL

#
# added restart, just because it's always a good idea
#
  restart;

#
# Added this just so that worksheet will display
# correctly ion the MApleprimes website
#
  interface(rtablesize=10):

de1 := (1-p)*(diff(f(x), `$`(x, 3)))+p*(diff(f(x), `$`(x, 3))+(1/2)*f(x)*(diff(f(x), `$`(x, 2))))

(1-p)*(diff(diff(diff(f(x), x), x), x))+p*(diff(diff(diff(f(x), x), x), x)+(1/2)*f(x)*(diff(diff(f(x), x), x)))

(1)

de2 := (1-p)*(diff(g(x), `$`(x, 2)))/Pr+p*((diff(g(x), `$`(x, 2)))/Pr+(1/2)*f(x)*(diff(g(x), x)))

(1-p)*(diff(diff(g(x), x), x))/Pr+p*((diff(diff(g(x), x), x))/Pr+(1/2)*f(x)*(diff(g(x), x)))

(2)

ibvc := f(0), (D(f))(0), (D(f))(5)-1, g(0)-1, g(5); n := 5; F := unapply(add(b[k](x)*p^k, k = 0 .. n), x); G := unapply(add(c[k](x)*p^k, k = 0 .. n), x)

f(0), (D(f))(0), (D(f))(5)-1, g(0)-1, g(5)

 

5

 

proc (x) options operator, arrow; b[0](x)+b[1](x)*p+b[2](x)*p^2+b[3](x)*p^3+b[4](x)*p^4+b[5](x)*p^5 end proc

 

proc (x) options operator, arrow; c[0](x)+c[1](x)*p+c[2](x)*p^2+c[3](x)*p^3+c[4](x)*p^4+c[5](x)*p^5 end proc

(3)

DE1 := series(eval(de1, f = F), p = 0, n+1)

series(diff(diff(diff(b[0](x), x), x), x)+(diff(diff(diff(b[1](x), x), x), x)+(1/2)*b[0](x)*(diff(diff(b[0](x), x), x)))*p+(diff(diff(diff(b[2](x), x), x), x)+(1/2)*b[0](x)*(diff(diff(b[1](x), x), x))+(1/2)*b[1](x)*(diff(diff(b[0](x), x), x)))*p^2+(diff(diff(diff(b[3](x), x), x), x)+(1/2)*b[0](x)*(diff(diff(b[2](x), x), x))+(1/2)*b[1](x)*(diff(diff(b[1](x), x), x))+(1/2)*b[2](x)*(diff(diff(b[0](x), x), x)))*p^3+(diff(diff(diff(b[4](x), x), x), x)+(1/2)*b[0](x)*(diff(diff(b[3](x), x), x))+(1/2)*b[1](x)*(diff(diff(b[2](x), x), x))+(1/2)*b[2](x)*(diff(diff(b[1](x), x), x))+(1/2)*b[3](x)*(diff(diff(b[0](x), x), x)))*p^4+(diff(diff(diff(b[5](x), x), x), x)+(1/2)*b[0](x)*(diff(diff(b[4](x), x), x))+(1/2)*b[1](x)*(diff(diff(b[3](x), x), x))+(1/2)*b[2](x)*(diff(diff(b[2](x), x), x))+(1/2)*b[3](x)*(diff(diff(b[1](x), x), x))+(1/2)*b[4](x)*(diff(diff(b[0](x), x), x)))*p^5+O(p^6),p,6)

(4)

DE2 := series(eval(de2, g = G), p = 0, n+1)

series((diff(diff(c[0](x), x), x))/Pr+((diff(diff(c[1](x), x), x))/Pr+(1/2)*f(x)*(diff(c[0](x), x)))*p+((diff(diff(c[2](x), x), x))/Pr+(1/2)*f(x)*(diff(c[1](x), x)))*p^2+((diff(diff(c[3](x), x), x))/Pr+(1/2)*f(x)*(diff(c[2](x), x)))*p^3+((diff(diff(c[4](x), x), x))/Pr+(1/2)*f(x)*(diff(c[3](x), x)))*p^4+((diff(diff(c[5](x), x), x))/Pr+(1/2)*f(x)*(diff(c[4](x), x)))*p^5+O(p^6),p,6)

(5)

CO := map(coeffs, eval([ibvc], f = F), p); CT := map(coeffs, eval([ibvc], g = G), p)

[b[0](0), b[1](0), b[2](0), b[3](0), b[4](0), b[5](0), (D(b[0]))(0), (D(b[1]))(0), (D(b[2]))(0), (D(b[3]))(0), (D(b[4]))(0), (D(b[5]))(0), (D(b[0]))(5)-1, (D(b[1]))(5), (D(b[2]))(5), (D(b[3]))(5), (D(b[4]))(5), (D(b[5]))(5), g(0)-1, g(5)]

 

[f(0), (D(f))(0), (D(f))(5)-1, c[0](0)-1, c[1](0), c[2](0), c[3](0), c[4](0), c[5](0), c[0](5), c[1](5), c[2](5), c[3](5), c[4](5), c[5](5)]

(6)

NULL

for k from 0 to n do IBVC1 := select(has, CO, b[k]); slv := dsolve({coeff(DE1, p, k), op(IBVC1)}); b[k] := unapply(rhs(slv), x) end do

F(x) = F(x)+O(p^(n+1))

(1/10)*x^2+(-(1/6000)*x^5+(5/96)*x^2)*p+((11/20160000)*x^8-(1/5760)*x^5+(325/16128)*x^2)*p^2+(-(1/532224000)*x^11+(11/12902400)*x^8-(29/258048)*x^5+(3125/1548288)*x^2)*p^3+((9299/1452971520000000)*x^14-(1/255467520)*x^11+(671/867041280)*x^8-(775/18579456)*x^5-(157028125/37196070912)*x^2)*p^4+(-(1272379/59281238016000000000)*x^17+(9299/557941063680000)*x^14-(157/34334834688)*x^11+(5665/12485394432)*x^8+(484625/127529385984)*x^5-(11930796875/3570822807552)*x^2)*p^5 = (1/10)*x^2+(-(1/6000)*x^5+(5/96)*x^2)*p+((11/20160000)*x^8-(1/5760)*x^5+(325/16128)*x^2)*p^2+(-(1/532224000)*x^11+(11/12902400)*x^8-(29/258048)*x^5+(3125/1548288)*x^2)*p^3+((9299/1452971520000000)*x^14-(1/255467520)*x^11+(671/867041280)*x^8-(775/18579456)*x^5-(157028125/37196070912)*x^2)*p^4+(-(1272379/59281238016000000000)*x^17+(9299/557941063680000)*x^14-(157/34334834688)*x^11+(5665/12485394432)*x^8+(484625/127529385984)*x^5-(11930796875/3570822807552)*x^2)*p^5+O(p^6)

(7)

plot(eval(F(x), p = 1), x = 0 .. 5, color = blue); plot(eval(diff(F(x), x), p = 1), x = 0 .. 5, color = red)

 

 

for k from 0 by .2 to 5 do NM1 := subs[eval](p = 1, F(k)) end do

``

``

for k from 0 to 1 do IBVC2 := select(has, CT, c[k]); sol := dsolve({coeff(DE2, p, k), op(IBVC2)}); c[k] := unapply(rhs(sol), x) end do

[c[0](0)-1, c[0](5)]

 

c[0](x) = -(1/5)*x+1

 

proc (x) options operator, arrow; -(1/5)*x+1 end proc

 

[c[1](0), c[1](5)]

 

c[1](x) = Int(Int((1/10)*f(_z1)*Pr, _z1 = 0 .. _z1), _z1 = 0 .. x)-(1/5)*(Int(Int((1/10)*f(_z1)*Pr, _z1 = 0 .. _z1), _z1 = 0 .. 5))*x

 

proc (x) options operator, arrow; Int(Int((1/10)*f(_z1)*Pr, _z1 = 0 .. _z1), _z1 = 0 .. x)-(1/5)*(Int(Int((1/10)*f(_z1)*Pr, _z1 = 0 .. _z1), _z1 = 0 .. 5))*x end proc

(8)

 

NULL

Download hpm2.mw

it is possible to use the 'parameters' option with the dsolve/numeric command.

This will produce a solution module with explicit values of 'p', 'lambda' and 'rho' as yet unspecified

These parameters can be subsequently given values, and relevant outputs obtained, as in the following.

Note that depending on the parameter values which are used, the ode system may (or may not) exhibit a singularity in the solution process

  restart;

  with(plots):
  interface(rtablesize=10):

  N:= 10:
  beta:= 1:
  alpha:= (N + 2)/(N - 2):
  ics:= {U(0) = beta, V(0) = 1, X(0) = 0, Y(0) = 1}:
  sys_ode:= { diff(U(r), r) = r^(N - 1)*rho^(N - 2)*p*X(r),
              diff(V(r), r) = r^(N - 1)*rho^(N - 2)*p*Y(r),
              diff(X(r), r) = -r^(N - 1)*rho^N*(U(r)^(alpha - 1) + lambda*U(r)),
              diff(Y(r), r) = -r^(N - 1)*rho^N*(lambda + (2^alpha - 1)*U(r)^(alpha - 2))*V(r)
            }:
  sol:= dsolve( [ sys_ode[], ics[] ],
                  numeric,
                  parameters = [lambda, p, rho]
              ):

#
# List the unknown parameters (just as a check!)
#
  sol(parameters);
#
# Supply one set of parameter values. Note
# this is dependent on order.
#
# Then plot the solution
#
  sol( parameters = [1, 1, 1] ):
  odeplot( sol,
           [ [r, U(r)], [r, V(r)], [r,X(r)], [r, Y(r)] ],
             r=0..1,
             color=[red, blue, green ,black]
         );
#
# Supply another set of parameter values and
# plot the solution
#
  sol( parameters = [2, 2, 2] ):
  odeplot( sol,
           [ [r, U(r)], [r, V(r)], [r,X(r)], [r, Y(r)]],
             r=0..1,
             color=[red, blue, green ,black]
         );
#
# Supply another set of parameter values and
# plot the solution
#
  sol( parameters = [1, 2, 3] ):
  odeplot( sol,
           [ [r, U(r)], [r, V(r)], [r,X(r)], [r, Y(r)]],
           r=0..1,
           color=[red, blue, green ,black]
         );
                  

[lambda = undefined, p = undefined, rho = undefined]

 

 

Warning, cannot evaluate the solution further right of .64415942, probably a singularity

 

 

Warning, cannot evaluate the solution further right of .45598040, probably a singularity

 

 

 

NULL


 

Download odePars.mw

.

2 3 4 5 6 7 8 Last Page 4 of 103