tomleslie

12497 Reputation

19 Badges

12 years, 231 days

MaplePrimes Activity


These are replies submitted by tomleslie

@planetmknzm 

rather than some "excerpt" which you *think* is appropriate. You can upload a worksheet using the big green up-arrow in the Mapleprimes toolbar. That way everybody here gets to see exactly what you are doing, rather than random rather meaningless code

For illustration of the problem whihc you are making impossible to diagnose, consider the first two lines of the "code" which you post, that is

c:=sphere([0,0,0],1,color=white,transparency=0.5):             
h3:=spacecurve([-x,-y,z],omega=ws..0,color=green,thickness=5):

Now let us consider the problems in these two line alone

  1. The sphere() command can only be executed after a with(plottools) command, so this command in your supplied command will never execute, it will just return as output sphere([0,0,0],1,color=white,transparency=0.5)
  2. The spacecurve() command can only be executed after a with(plots) command, so this command in your supplied command will never execute it will just return as output spacecurve([-x,-y,z],omega=ws..0,color=green,thickness=5):
  3. Even after loading the 'plots' package, the above spacecurve command will not execute because
    1. The value of 'ws' is undefined - it needs to evaluate to a numeric quantity which you do not supply
    2. The spacecurve parameter 'omega' is not used in the calculation of the spacecurve coordinates - unless you have a definition of the expressions 'x', 'y', 'z' in terms of omega - which you do not supply
  4. After two lines of so-called code, your worksheet is useless - I haven't even bthered to look at the rest of it
  5. Upload the code you are actually running

 

@Saha 

which is posted again here for your convenience. The only difference between the attached and my original is that I have included the command interface(version) immediately after the restart, just to demonstrate that this executes correctly in your version of Maple.

I fixed all your problems once. I can't think of any reason why I should do the same again

  restart;
  interface(version);
  with(plots):

`Standard Worksheet Interface, Maple 18.02, Windows 7, October 20 2014 Build ID 991181`

(1)

  eq1:= diff(f(x), x$3)
        +
        1/2*(1-phi)^2.5*(1-phi+phi*rho__s/rho__f1)*(eta*cos__omega + f(x)*sin__omega)*diff(f(x), x$2)
        +
        (1-phi)^2.5*M*sin(alpha)^2*(1-diff(f(x),x))
        +
        (1-phi)^2.5*(1-phi+phi*rho__beta__s/rho__beta__f1)*lambda__T*theta(x);

  eq2:= K__nf/K__f*diff(theta(x),x$2)
        +
        Pr/2*(eta*cos__omega + f(x)*sin__omega)*diff(theta(x),x);
#
# The OP's original list of boundary conditions contained f(1)=0.
# In a 1-dimensional BVP only two boundaries are allowed, so the
# condition f(1)=0 has been "replaced" by the condition
#
#    D(f)(0) = val
#
# A "shooting" method will be used later to determine a value for
# "val", whihc ensures that the original BC ( ie f(1)=0 ) is
# satisfied
#
  bcs := f(0) = 0, D(f)(0) = val, f(10) = 1, theta(0) = 1, theta(10) = 0;

  params:= [ K__f=0.613, K__nf=0.6842, M=1, Pr=6.2,
             alpha=0, eta=0, phi=0.05, rho__f1=997.1,
             rho__s=5200, cos__omega=1, lambda__T=0,
             rho__beta__f1=20939.1, rho__beta__s=997.1,
             sin__omega=1
           ];
#
# So the ODE system with all parameter values inserted is
#
  odesys:=eval([eq1, eq2], params);

diff(diff(diff(f(x), x), x), x)+(1/2)*(1-phi)^2.5*(1-phi+phi*rho__s/rho__f1)*(eta*cos__omega+f(x)*sin__omega)*(diff(diff(f(x), x), x))+(1-phi)^2.5*M*sin(alpha)^2*(1-(diff(f(x), x)))+(1-phi)^2.5*(1-phi+phi*rho__beta__s/rho__beta__f1)*lambda__T*theta(x)

 

K__nf*(diff(diff(theta(x), x), x))/K__f+(1/2)*Pr*(eta*cos__omega+f(x)*sin__omega)*(diff(theta(x), x))

 

f(0) = 0, (D(f))(0) = val, f(10) = 1, theta(0) = 1, theta(10) = 0

 

[K__f = .613, K__nf = .6842, M = 1, Pr = 6.2, alpha = 0, eta = 0, phi = 0.5e-1, rho__f1 = 997.1, rho__s = 5200, cos__omega = 1, lambda__T = 0, rho__beta__f1 = 20939.1, rho__beta__s = 997.1, sin__omega = 1]

 

[diff(diff(diff(f(x), x), x), x)+.5325197465*f(x)*(diff(diff(f(x), x), x)), 1.116150082*(diff(diff(theta(x), x), x))+3.100000000*f(x)*(diff(theta(x), x))]

(2)

#
# The shooting method
#
# This determines a value for D(f)(0) such that the
# OP's original desired condition f(1)=0 is satisfied
#
  SM:= proc( shoot )
             local ans:
             if type( shoot, numeric)
             then ans:= dsolve
                        ( eval
                          ( [odesys[], bcs],
                            val=shoot
                          ),
                          numeric
                        ):
                  return rhs(ans(1)[2]);
             else return 'procname(_passed)'
             fi;
       end proc:
  sv:= fsolve( SM );

-0.1290531852e-1

(3)

#
# Solve the system
#
  sol:= dsolve
        (  eval
           ( [odesys[], bcs],
             val=sv
           ),
          numeric
        ):
#
# Check that when x=1, f(x)=0, which was OP's desired condition
#
  sol(1)[1..2];
#
# Plot the results
#
  odeplot
  ( sol,
    [ [x, f(x)],
      [x, theta(x)]
    ],
    x = 0 .. 10,
    color=[red, blue],
    legend=[typeset(f(x)), typeset(theta(x))]
  );
  

[x = 1., f(x) = HFloat(1.9274077505456594e-12)]

 

 

 


 

Download odeProb.mw

 

@Tamour_Zubair

With N=5 and t=0..3, your dependent varaibles have numerical values in the range +/-10^27.

With N=6 and t=0..3, your dependent varaibles have numerical values in the range +/-10^45

With N=7 and t=0..3, you dependent varaibles have numerical values in the range +/-10^63

So every time you increase the number of equations in your system, the numerical range of the dependend variables increases by about 10^18. Based on the above, with N=8, one might expect the dependent variables to have numerical values in the range +/-10^81, which is about the same as the number of atoms in the known universe

I always get worried when I see numbers like this. Ask yourself some serious questions.

  1. If you are modelling some "real-world" phenomenon, what range of values would you expect for the dependent variables? You ought to be able to make a rough guess!
  2. Do you really believe that your system can produce values greater than the number of atoms in the universe? If so, how?, why?, what are you modelling?

@Tamour_Zubair 

You state

Thanks a lot for this. I cannot use builton command beacuse system of ODE's increase as i will increase the value of "N".

For a start, there is no 'N' in your worksheet, so what is increasing? The number of ODEs? the complexity of the ODEs?

Why do you think that you can produce some "magic code" which will work better than Maple's built-in ODE solvers?

Maybe if I understood the real problem you are trying to solve, I could make sensible suggestions

@planetmknzm 

using the big green up-arrow in the Mapleprimes toolbar to upload worksheets which are causing you a problem?

That way, responders here can execute your worksheet, figure out  possible improvements, and do all kinds of useful stuff.

On the other hand you can upload a video which doesn't show (most of) your code, so cannot be read , cannot be executed, cannot be analysed, and hence is unlikely to be "improved"

  1. Rule 1 - upload an executable worksheet
  2. Rule 2 - if you think there is a "better" approach, stop thinking and apply Rule 1

@Kitonum 

You could set the animation speed using the FPS box in Maple's animation toolbar

you still don't need the InertForm() package, entering data as strings, then parsing strings etc etc - Sure it will work, but that is the conplicated way to do something simple. See the attached for a "simple" variation on my original answer, based on the following observations

  1. entering data as strings - cumbersome and unnecessary
  2. "parsing" such entries - cumbersome and unnecessary
  3. of course the above approaches will "work", but why bother? Isn't the attached just - well, simpler?

#
# One possible variant of OP's desired expression
#
  Set:= `%+`(  1/3%*<3^(1/2), -3^(1/2), 3^(1/2), 0>,
              -1/3%*<3^(1/2),  3^(1/2), 3^(1/2), 0>
              +3*<2,1,1,0>
            );
#
# Quick check on export+import to/from MathML
#
  a:= MathML:-Export(Set):
  MathML:-Import(a);

`%+`(`%*`(1/3, Vector(4, {(1) = sqrt(3), (2) = -sqrt(3), (3) = sqrt(3), (4) = 0})), -`%*`(1/3, Vector(4, {(1) = sqrt(3), (2) = sqrt(3), (3) = sqrt(3), (4) = 0}))+(Vector(4, {(1) = 6, (2) = 3, (3) = 3, (4) = 0})))

 

`%+`(`%*`(1/3, Vector[column](%id = 36893488148080653604)), -`%*`(1/3, Vector[column](%id = 36893488148080653724))+Vector[column](%id = 36893488148080653844))

(1)

#
# Another possible variant of OP's desired expression
#
  Set:= `%+`(  1/3%*<3^(1/2), -3^(1/2), 3^(1/2), 0>,
              -1/3%*<3^(1/2),  3^(1/2), 3^(1/2), 0>,
               <2,1,1,0>
            );
#
# Quick check on  export+import to/from MathML
#
  a:= MathML:-Export(Set):
  MathML:-Import(a);

`%+`(`%*`(1/3, Vector(4, {(1) = sqrt(3), (2) = -sqrt(3), (3) = sqrt(3), (4) = 0})), -`%*`(1/3, Vector(4, {(1) = sqrt(3), (2) = sqrt(3), (3) = sqrt(3), (4) = 0})), Vector(4, {(1) = 2, (2) = 1, (3) = 1, (4) = 0}))

 

`%+`(`%*`(1/3, Vector[column](%id = 36893488148080646372)), -`%*`(1/3, Vector[column](%id = 36893488148080646492)), Vector[column](%id = 36893488148080646612))

(2)

 

Download inert3.mw

 

whether you use "worksheet mode", "Document mode", "1-D Input", "2D-input" what. You can save the worksheet on your local machine as a file with extension ".mw". Upload this file to this site using the big green up-arrow in the Mapleprimes toolbar

post an executable worksheet using the big green up-arrow in the Mapleprimes toolbar - because no-one here is going to re-type your stuff from a useless "picture"

take a long hard look at the boundary condition

theta(0, t) = 1 + b(1 - cos(varpi*t))

'b' has been defined as 1, so the above is equivalent to

theta(0, t) = 1 + 1(1 - cos(varpi*t))

The second term on the rhs will be interpreted as as a function named '1' with an argument (which is irrelevant) because in Maple the construction 1(x) will return 1, so the above bc evaluates to

theta(0, t) = 2

I'm guessing that this isn't what you intended!

@Rouben Rostamian  

which applies to the response you have given, as well as that gien by dharr.

(And before I make this observation, I want to make it clear that I have been playing with this problem for a while, and so far haven't come with anything significant)

My issue with both Rouben's response (and that of dharr), is that both imply that the velocity in the x-direction changes. My simple-mided brain says that there is no "force" in the x-direction, hence there can be no acceleration in the x-direction, hence the x-velocity cannot change.

So the question for boith Rouben and dharr - why does the x-velocity change? What force is causing this?

as shown in the attached

  restart;
  expr:= diff(u(t,x,v),t)+v*diff(u(t,x,v),x)+diff(u(t,x,v),x)*int(u(t,x,v)*v*exp(x),v)-v*diff(u(t,x,v),v)*sin(v)*int(v*diff(u(t,x,v),x),v)=0;
  kf:=convert(FunctionAdvisor(known_functions, quiet), set):
  seq( `if`(member(op(0,j), kf), j, NULL), j in indets(expr, function));

diff(u(t, x, v), t)+v*(diff(u(t, x, v), x))+(diff(u(t, x, v), x))*(int(u(t, x, v)*v*exp(x), v))-v*(diff(u(t, x, v), v))*sin(v)*(int(v*(diff(u(t, x, v), x)), v)) = 0

 

exp(x), sin(v)

(1)

 

Download isol2.mw

@Luxio 

my first answer still applies!

However the following "toy" example shows usage of both the CodeGeneration:-Python() command and the Python() package, just to give some idea of what is possible.

  restart:

  alias(CG=CodeGeneration):
  alias(P=Python):
#
# Define a not-too-trivial" procedure in Maple
#
  detHilbert:=  proc(n::posint)
                    uses LinearAlgebra;
                    return Determinant( HilbertMatrix( n ) );
                end proc:
#
# Compute one value for later comparison purposes
#
  evalf(detHilbert(10));

0.2164179226e-52

(1)

#
# Generate "Python code" for the above procedure (as a string)
#
  pycode:= CodeGeneration:-Python(detHilbert, output=string);

"import numpy.linalg
import scipy.linalg

def detHilbert (n):
    return(numpy.linalg.det(scipy.linalg.hilbert(n)))
"

(2)

#
# Define the above procedure within Python
#
  Python:-EvalString( pycode, output=none );
#
# Execute the Python procedure (in Python)
#
  Python:-EvalString("detHilbert(10)")

HFloat(2.1642457579118234e-53)

(3)

 

Download pyMaple.mw

Maple 2022.1 produces the attached - which at least avoids the Float(undefined) issue


 

NULL

restart

interface(version)

`Standard Worksheet Interface, Maple 2022.1, Windows 7, May 26 2022 Build ID 1619613`

(1)

with(Student[NumericalAnalysis])

NULL

fff := (`&eta;__1`^2-1.)^2*(`&zeta;__1`^2-1.)^2*(0.12053033309e-3*`&zeta;__1`-4.5961086182*10^(-8)*`&eta;__1`^2*`&zeta;__1`^4-5.6584306313*10^(-8)*`&eta;__1`^4*`&zeta;__1`^4-1.5606093643*10^(-7)*`&zeta;__1`*`&eta;__1`^2-4.168658171*10^(-7)*`&zeta;__1`*`&eta;__1`-1.6404108106*10^(-7)*`&zeta;__1`^2*`&eta;__1`-1.0454641424*10^(-8)*`&eta;__1`*`&zeta;__1`^4+1.0615594865*10^(-8)*`&eta;__1`^5*`&zeta;__1`^4+5.4851856568*10^(-9)*`&eta;__1`^3*`&zeta;__1`^2+5.6161016651*10^(-9)*`&eta;__1`^3*`&zeta;__1`^4-1.4132765167*10^(-8)*`&eta;__1`^5*`&zeta;__1`^2-7.8157365683*10^(-7)*`&zeta;__1`*`&eta;__1`^4+2.9373057668*10^(-8)*`&eta;__1`^5*`&zeta;__1`^5-5.2577413654*10^(-8)*`&eta;__1`^5*`&zeta;__1`^3-7.032574429*10^(-8)*`&eta;__1`^3*`&zeta;__1`^5+2.3272955826*10^(-8)*`&eta;__1`^5*`&zeta;__1`+1.5782217112*10^(-7)*`&eta;__1`^3*`&zeta;__1`^3+2.1771522925*10^(-8)*`&eta;__1`*`&zeta;__1`^5-8.3051507888*10^(-8)*`&eta;__1`^3*`&zeta;__1`-2.2997952126*10^(-8)*`&eta;__1`*`&zeta;__1`^3+0.22608138853e-5*`&eta;__1`^4-0.10412827312e-5*`&zeta;__1`^4-0.53422910417e-5*`&zeta;__1`^3-0.53519692056e-5*`&eta;__1`^2+0.54946587424e-4*`&zeta;__1`^2+7.5007400675*10^(-7)*`&zeta;__1`^5-6.6427826628*10^(-9)*`&eta;__1`^3+3.5821686059*10^(-9)*`&eta;__1`^5-2.4361928132*10^(-7)*`&eta;__1`+0.84495460425e-5*`&eta;__1`^2*`&zeta;__1`^2-0.44664222239e-5*`&eta;__1`^4*`&zeta;__1`^2-0.4071733639e-5*`&eta;__1`^4*`&zeta;__1`^3+0.87015148204e-5*`&eta;__1`^2*`&zeta;__1`^3+0.25592467903e-5*`&eta;__1`^4*`&zeta;__1`^5-0.5430102455e-5*`&eta;__1`^2*`&zeta;__1`^5+0.62041471332e-4)/(-0.22768667085e-2*`&eta;__1`-.20653118433*`&zeta;__1`+0.11770764739e-2*`&eta;__1`^9*`&zeta;__1`^2+.14111119517*`&eta;__1`^4+.66506159208*`&zeta;__1`^4-.60151130424*`&zeta;__1`^5+.58641863992*`&zeta;__1`^3+0.49423587807e-2*`&eta;__1`^3-0.32350168299e-2*`&eta;__1`^5-1.082755589*`&eta;__1`^2-1.3163499567*`&zeta;__1`^2+0.52300044044e-1*`&eta;__1`^5*`&zeta;__1`^7-0.40389225121e-1*`&eta;__1`^3*`&zeta;__1`^7+0.10746052526e-1*`&eta;__1`*`&zeta;__1`^7+0.37597513815e-2*`&eta;__1`^9*`&zeta;__1`^7+0.26109384594e-2*`&eta;__1`^7*`&zeta;__1`^9-0.70972300998e-2*`&eta;__1`^5*`&zeta;__1`^9+0.63616448215e-2*`&eta;__1`^3*`&zeta;__1`^9-0.18753531811e-2*`&eta;__1`*`&zeta;__1`^9-.22701198791*`&eta;__1`^6*`&zeta;__1`^6+0.56271163518e-2*`&eta;__1`^4*`&zeta;__1`^8-0.48478855017e-1*`&eta;__1`^8*`&zeta;__1`^4-.22761063366*`&zeta;__1`^6*`&eta;__1`^2+.46301464814*`&zeta;__1`^6*`&eta;__1`^4-0.18114547653e-2*`&eta;__1`^9*`&zeta;__1`^4+0.81527768559e-3*`&eta;__1`^9*`&zeta;__1`^6+.65574325943-.15635457174*`&eta;__1`^8*`&zeta;__1`+.64029273264*`&eta;__1`^8*`&zeta;__1`^3-1.8403657443*`&eta;__1`^6*`&zeta;__1`^7-0.22007436935e-2*`&eta;__1`^2*`&zeta;__1`^8-0.30978282387e-2*`&zeta;__1`*`&eta;__1`+.70551660939*`&zeta;__1`*`&eta;__1`^2+0.38970885732e-2*`&zeta;__1`^2*`&eta;__1`-.62093209055*`&eta;__1`^2*`&zeta;__1`^4+1.9334990569*`&eta;__1`^2*`&zeta;__1`^2-.80179945016*`&eta;__1`^4*`&zeta;__1`^4+.1920464905*`&eta;__1`^4*`&zeta;__1`^2-0.81381430978e-3*`&eta;__1`*`&zeta;__1`^4+0.46097207851e-2*`&eta;__1`^3*`&zeta;__1`^4-0.93963570584e-2*`&eta;__1`^3*`&zeta;__1`^2-0.85894534062e-2*`&eta;__1`^5*`&zeta;__1`^4+0.8278524871e-2*`&eta;__1`^5*`&zeta;__1`^2+4.2209528188*`&eta;__1`^4*`&zeta;__1`^3-2.669958003*`&eta;__1`^2*`&zeta;__1`^3-.94779423754*`&zeta;__1`*`&eta;__1`^4-5.9197768267*`&eta;__1`^4*`&zeta;__1`^5+3.4563944947*`&eta;__1`^2*`&zeta;__1`^5-.10468335582*`&eta;__1`^5*`&zeta;__1`^5+0.80855499912e-1*`&eta;__1`^5*`&zeta;__1`^3+0.7601825081e-1*`&eta;__1`^3*`&zeta;__1`^5-0.21374958034e-1*`&eta;__1`^5*`&zeta;__1`-0.56315405543e-1*`&eta;__1`^3*`&zeta;__1`^3-0.18963873748e-1*`&eta;__1`*`&zeta;__1`^5+0.14324735033e-1*`&eta;__1`^3*`&zeta;__1`+0.13191002642e-1*`&eta;__1`*`&zeta;__1`^3-0.301782967e-2*`&eta;__1`^6*`&zeta;__1`^8+.34307133883*`&eta;__1`^6+.26989142602*`&zeta;__1`^7+0.75042415188e-3*`&eta;__1`^7-0.57170204466e-1*`&eta;__1`^8-0.40854298828e-3*`&zeta;__1`^8-0.18089939414e-3*`&eta;__1`^9-0.48267577378e-1*`&zeta;__1`^9-0.40463518464e-2*`&zeta;__1`^6-0.26416622831e-1*`&eta;__1`^7*`&zeta;__1`^7-0.20189726843e-2*`&eta;__1`^9*`&zeta;__1`-.81152175007*`&eta;__1`^8*`&zeta;__1`^5+.32758358916*`&eta;__1`^8*`&zeta;__1`^7+.60516338422*`&eta;__1`^6*`&zeta;__1`-2.7777061884*`&eta;__1`^6*`&zeta;__1`^3+3.8764153863*`&eta;__1`^6*`&zeta;__1`^5+0.64408301026e-3*`&eta;__1`^3*`&zeta;__1`^8+0.14976271107e-3*`&eta;__1`*`&zeta;__1`^8-0.17374541537e-2*`&eta;__1`^5*`&zeta;__1`^8+0.77976967502e-2*`&eta;__1`^9*`&zeta;__1`^3-0.95384754473e-2*`&eta;__1`^9*`&zeta;__1`^5+0.9436084324e-3*`&eta;__1`^7*`&zeta;__1`^8-0.95617026598e-3*`&eta;__1`*`&zeta;__1`^6-0.79980551762e-3*`&eta;__1`^3*`&zeta;__1`^6+0.52833995188e-2*`&eta;__1`^5*`&zeta;__1`^6-0.43427014208e-2*`&eta;__1`^7*`&zeta;__1`^6+0.66050016962e-2*`&eta;__1`^7*`&zeta;__1`^4-0.39563328597e-2*`&eta;__1`^7*`&zeta;__1`^2-.9191903249*`&eta;__1`^6*`&zeta;__1`^2+.80614880365*`&eta;__1`^6*`&zeta;__1`^4+2.9678721471*`&zeta;__1`^7*`&eta;__1`^4-1.724981418*`&zeta;__1`^7*`&eta;__1`^2+.13649316215*`&eta;__1`^6*`&zeta;__1`^9+.23302831691*`&eta;__1`^2*`&zeta;__1`^9-.32125390168*`&eta;__1`^4*`&zeta;__1`^9+0.12167023924e-1*`&eta;__1`^7*`&zeta;__1`-0.45528793761e-1*`&eta;__1`^7*`&zeta;__1`^3+0.57167454208e-1*`&eta;__1`^7*`&zeta;__1`^5+.10999473421*`&eta;__1`^8*`&zeta;__1`^2-0.43456747248e-2*`&eta;__1`^8*`&zeta;__1`^6)

plot3d(sqrt(fff), `&zeta;__1` = -1 .. 1, `&eta;__1` = -1 .. 1, color = green)

 

Quadrature(Quadrature(sqrt(fff), `&zeta;__1` = -1 .. 1, method = romberg[3]), `&eta;__1` = -1 .. 1, method = romberg[3])

0.2745463666e-1+0.*I

(2)

Quadrature(Quadrature(sqrt(fff), `&zeta;__1` = -1 .. 1, method = romberg[4]), `&eta;__1` = -1 .. 1, method = romberg[4])

0.3314502549e-1+0.*I

(3)

Quadrature(Quadrature(sqrt(fff), `&zeta;__1` = -1 .. 1, method = romberg[5]), `&eta;__1` = -1 .. 1, method = romberg[5])

0.3621732017e-1+0.*I

(4)

Quadrature(Quadrature(sqrt(fff), `&zeta;__1` = -1 .. 1, method = romberg[6]), `&eta;__1` = -1 .. 1, method = romberg[6])

0.3781929307e-1+0.1394002550e-4*I

(5)

``


 

Download undef.mw

 

@rashmi 

Even for j=1, lets say red colour curve, there are two curves there

As I explained before, technically there is one curve, with a gap in it because the relevant equation for yj(r) generates complex values.

To see onlyone section of each curve , you can plot with r=0..10 or r=-10..0, rather than r=-10..10 as shown in the attached.

  restart;
  expr:=[ [ (1.428571429*(r^2+.49-4*r*(r^2+.49-2*exp(-.1/r))/(2*r-.2*exp(-.1/r)/r^2)))*csc((1/3)*Pi),
          sqrt(16*r^2*(r^2+.49-2*exp(-.1/r))/(2*r-.2*exp(-.1/r)/r^2)^2+.49*cos((1/3)*Pi)^2-2.040816327*(r^2+.49-4*r*(r^2+.49-2*exp(-.1/r))/(2*r-.2*exp(-.1/r)/r^2))^2*cot((1/3)*Pi)^2)
        ],
        [(1.428571429*(r^2+.49-4*r*(r^2+.49-2*exp(-.1/r))/(2*r-.2*exp(-.1/r)/r^2)))*csc((1/3)*Pi),
          -sqrt(16*r^2*(r^2+.49-2*exp(-.1/r))/(2*r-.2*exp(-.1/r)/r^2)^2+.49*cos((1/3)*Pi)^2-2.040816327*(r^2+.49-4*r*(r^2+.49-2*exp(-.1/r))/(2*r-.2*exp(-.1/r)/r^2))^2*cot((1/3)*Pi)^2)
        ],
        [(5.*(r^2+0.4e-1-4*r*(r^2+0.4e-1-2*exp(-.2/r))/(2*r-.2*exp(-.2/r)/r^2)))*csc((1/3)*Pi),
          sqrt(16*r^2*(r^2+0.4e-1-2*exp(-.2/r))/(2*r-.2*exp(-.2/r)/r^2)^2+0.1e-1*cos((1/3)*Pi)^2-25.*(r^2+0.4e-1-4*r*(r^2+0.4e-1-2*exp(-.2/r))/(2*r-.2*exp(-.2/r)/r^2))^2*cot((1/3)*Pi)^2)
        ],
        [(5.*(r^2+0.4e-1-4*r*(r^2+0.4e-1-2*exp(-.2/r))/(2*r-.2*exp(-.2/r)/r^2)))*csc((1/3)*Pi),
         -sqrt(16*r^2*(r^2+0.4e-1-2*exp(-.2/r))/(2*r-.2*exp(-.2/r)/r^2)^2+0.1e-1*cos((1/3)*Pi)^2-25.*(r^2+0.4e-1-4*r*(r^2+0.4e-1-2*exp(-.2/r))/(2*r-.2*exp(-.2/r)/r^2))^2*cot((1/3)*Pi)^2)
        ]
       ]:

 

#
# From the above, the first and second parametric
# curves will the have the same "x-value" and
#  "y-values" of opposite sign, so the second
# parametric curve is identical to the first
# parametric curve exceptfor a reflection in the x-axis.
#
# The same applies to the third and fourth parametric
# curves - identical except for reflection in x-axis
#
  plot( [seq( [expr[j][], r=-10..10], j=1..4)],
          -13 .. 5,
          -5 .. 5,
          color=[red, blue, green, black]);
  plot( [seq( [expr[j][], r=0..10], j=1..4)],
          -13 .. 5,
          -5 .. 5,
          color=[red, blue, green, black]);
  plot( [seq( [expr[j][], r=-10..0], j=1..4)],
         -13 .. 5,
         -5 .. 5,
         color=[red, blue, green, black]);

 

 

 

 

 

 

 

Download parPlot3.mw

 

3 4 5 6 7 8 9 Last Page 5 of 196