Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 181 days
Contact tomleslie

MaplePrimes Activity


These are replies submitted by tomleslie

As always...

tomleslie 13876
January 09 2019
0 0

it depends!

Somehow you have to "encapsulate" the code you want to execute multiple times. The two obvious ways to do this are

  1. A seq() or 'for' loop
  2. Enclose the original code in a proc()..end proc wrapper and call this procedure with loop index values.

Why don't you post your worksheet (using the big green up-arrow in the toolbar)? That way, someone here can probably figure out your best option

All communication...

tomleslie 13876
January 08 2019
0 0

on Maple problems takes place via this site

Do not send anything to my email address

Any problem you have with my previous answer, report it here

 

 

two reasons...

tomleslie 13876
January 06 2019
0 0

why this code will not run

  1. the procedure Anova2s uses a function called FCDF() which is nowhere defined
  2. the final  lprint() and printf() statements are completely garbled - unbalanced parentheses, inappropriate use of single quotes, etc, etc. Where did this code come from?

You could try...

tomleslie 13876
January 02 2019
0 0

@mehran rajabi 

reading my original response - carefully!!!

The system you define as case1 has

  1. Three ODEs in three dependent variables r[1](t), V(t), q(t) - This is good
  2. System is
    1. Degree 2 in r[1](t)  -> 2 ICs
    2. Degree 2 in q(t)  ->  2 ICs
    3. Degree 1 in V(t)  -> 1 IC

for a required total of five ICs - and you have five. This is good

Since you have the correct numbers of equations, dependent variables, and initial conditions, this system can (in principle) be solved. See my original post for the solution

###########################################################

The system you define as case2 has

  1. Three ODEs in four dependent variables r[1](t), r[2](t), V(t), q(t) - This is bad
  2. System is
    1. Degree 2 in r[1](t) -> 2 ICs
    2. Degree 2 in r[2](t) -> 2 ICs
    3. Degree 2 in q(t) ->  2 ICs
    4. Degree 1 in V(t) -> 1 IC

for a required total of seven ICs - but you only have five ICs. This is bad

Since you have the wrong numbers of equations, dependent variables, and initial conditions, this system cannot be solved.

Opops - Maple 18...

tomleslie 13876
December 30 2018
0 0

To work in Maple 18, my previous post requires a minor revision. See attached

> 

  restart;
  kernelopts(version);

`Maple 18.02, X86 64 WINDOWS, Oct 20 2014, Build ID 991181`

 (1)
> 

#
# Doing everything symbolically. NB if my interpretation
# of BCs, ICs is correct, then solution only exists in the
# range x=0..2
#
  restart;
  pde:=diff(u(x,t),t,t)=9*diff(u(x,t),x,x):
  bcs:=u(0,t)=0, u(2,t)=0:
  ics:= u(x,0) = piecewise( x<0,  undefined,
                            x<1,  x*(2-x),
                            x<=2, x,
                            undefined
                          ),
        D[2](u)(x,0) = 0:
  sol:=pdsolve( [pde, bcs, ics]);
  sol:=subs( _Z1=n, sol); # added this for Maple 18

u(x, t) = Sum(-4*(Pi^2*_Z1^2*(-1)^_Z1+Pi*sin((1/2)*Pi*_Z1)*_Z1+4*cos((1/2)*Pi*_Z1)-4)*sin((1/2)*Pi*_Z1*x)*cos((3/2)*Pi*_Z1*t)/(Pi^3*_Z1^3), _Z1 = 1 .. infinity)

  

u(x, t) = Sum(-4*(Pi^2*n^2*(-1)^n+Pi*sin((1/2)*Pi*n)*n+4*cos((1/2)*Pi*n)-4)*sin((1/2)*Pi*n*x)*cos((3/2)*Pi*n*t)/(Pi^3*n^3), n = 1 .. infinity)

 (2)
> 

#
# Since the soluton is an infinite sum with no
# closed form, this solution is "difficult" to
# plot. However one can easily extract the summand,
# and easily plot the first few terms of the sum
# ie n=1, n=2, n=3, n=4....
#
  solfn:=unapply( op([2,1],sol), [n,x,t]);
  nterms:= 4:
  for j from 1 by 1 to nterms do
      plot3d( solfn(j, x, t),
              x = 0..2,
              t = 0..10,
              title = typeset("n=", j," term"),
              titlefont = [times, bold, 20]
            )
  od;
#
# Alternatively, one could produce a plot of the sum
# of the first few terms
#
  nterms:= 50;
  plot3d( add
          ( solfn(j,x,t),
            j=1..nterms
          ),
          x = 0..2,
          t = 0..10,
          title = typeset("Sum of first ", nterms," terms"),
          titlefont = [times, bold, 20]
        );
#
# Or an animation, in time
#
  plots:-animate( plot,
           [ add
             ( solfn(j,x,t),
               j=1..nterms
             ),
             x = 0..2
           ],
           t = 0..10,
           frames=100
        );

proc (n, x, t) options operator, arrow; -4*(Pi^2*n^2*(-1)^n+Pi*sin((1/2)*n*Pi)*n+4*cos((1/2)*n*Pi)-4)*sin((1/2)*Pi*n*x)*cos((3/2)*Pi*n*t)/(Pi^3*n^3) end proc

  

 

 

 

 

50

  

 

 
> 

##########################################
# Doing everything numerically
#
  restart;
  np:=1000:
  pde:=diff(u(x,t),t,t)=9*diff(u(x,t),x,x):
  bcs:=u(0,t)=0, u(2,t)=0:
  ics:= u(x,0) = piecewise( x<1,  x*(2-x),
                            x<=2, x
                          ),
        D[2](u)(x,0) = 0:
  
  sol:= pdsolve( pde, [bcs, ics], numeric, time=t, range=0..2, timestep=0.01, spacestep=5/np);
  sol:-plot3d( t=0..10, x=0..2, numpoints=10000);
#
# Alternatively one can produce an animation in time
# using something like
#
  sol:-animate( t=10,
                frames=100,
                title="time = %f"
              );

module () local INFO; export plot, plot3d, animate, value, settings; option `Copyright (c) 2001 by Waterloo Maple Inc. All rights reserved.`; end module

  

 

 
> 

 


 

Download stringPDE3.mw

Errrr??...

tomleslie 13876
December 24 2018
0 0

dsolve() takes a 'method=something' option, but I wasn't aware that this option is available for pdsolve() - and I can't find any reference in the help that it does??? Can you provide reference to relevant pdsolve() help page

Yet more formatted print statements...

tomleslie 13876
December 23 2018
0 0

@cinderella 

See the attached, which does your latest request, two ways.

(Output format is neater in a worksheet than it appears on this site - honest)

> 

#
# Initialize and set up necessary functions
#
  restart;
  with(Optimization):
  y := (p, e) -> alpha-beta*p+k*e:
  g := e-> (1/2)*mu*e^2:
  Lambda:= z-> int((z-u)*phi(u), u=0..z):
  PI := unapply( (p-c)*(y(p, e)-z)-(p-v)*Lambda(z)-g(e), [p,e,z]):
  Pi_rd:= unapply( (p-w)*(y(p, e)-z)-(p-v)*Lambda(z)-g(e), [p, e, z]):
> 

#
# Set up k-values and optimizer conditions
#
  kVals:=[1, 2, 3]:
  optConds:=  p = 0 .. 60,
              e = -1000 .. 30,
              z = -10000 .. 25,
              initialpoint = {e = 0, p = 10, z = 0},
              maximize:

  for ind from 1 by 1 to numelems(kVals) do
    #
    # Set up parameters
    #
      pars:= [ k = kVals[ind],
               beta = 2,
               alpha = 50,
               mu = 10,
               c = 5,
               v = 1,
               phi(u) = 1/2,
               w=10
             ]:
      Int_Profit[ind]:= NLPSolve
                        ( eval
                          ( PI(p, e, z),
                            pars
                          ),
                          optConds
                       );
      qopt[ind]:= eval(z+y(p, e), [ Int_Profit[ind][2][], pars[]]);
      Re_Profit[ind]:= NLPSolve
                       ( eval
                         ( Pi_rd(p, e, z),
                           pars
                         ),
                         optConds
                       );
      qd_opt[ind]:= eval(z+y(p, e), [ Re_Profit[ind][2][], pars[]]);
      M_Profit[ind]:= eval((w-c)*qd_opt[ind], pars);
      Whole_Profit[ind]:= M_Profit[ind]+Re_Profit[ind][1];
      t[ind]:= eval((Int_Profit[ind][1]-Re_Profit[ind][1])/(w-c), pars);
  od:
> 

#
# How to access results, eg for Re_Profit.
#
# Others ie Int_Profit
#           Whole_Profit
#           qopt
#           qd_opt
#           t
#
# are accessed similarly
#
  Re_Profit[1];
  Re_Profit[2];
  Re_Profit[3];

[118.994532457413399, [e = HFloat(0.78757902976968), p = HFloat(17.875790315337227), z = HFloat(-0.9333832842273811)]]

  

[129.081143396166169, [e = HFloat(1.7075918734530875), p = HFloat(18.537959441958364), z = HFloat(-0.9736548233351248)]]

  

[150.271340626426877, [e = HFloat(2.9781181923084055), p = HFloat(19.9270613673776), z = HFloat(-1.0489806678261728)]]

 (1)
> 

  printf( "\n\n\t\t                 k=1           k=2           k=3\n\n");
  printf( "\t\t      qd_opt%14.8f%14.8f%14.8f\n",
          seq( qd_opt[j],
               j=1..numelems(kVals)
             )
        ):
  printf( "\t\t   Re_Profit%14.8f%14.8f%14.8f\n",
          seq( Re_Profit[j][1],
               j=1..numelems(kVals)
             )
        ):
  printf( "\t\t  Int_Profit%14.8f%14.8f%14.8f\n",
          seq( Int_Profit[j][1],
               j=1..numelems(kVals)
             )
        ):
  printf( "\t\tWhole_Profit%14.8f%14.8f%14.8f\n",
          seq( Whole_Profit[j],
               j=1..numelems(kVals)
             )
        ):
  printf( "\t\t           t%14.8f%14.8f%14.8f\n\n\n",
          seq( t[j],
               j=1..numelems(kVals)
             )
        );
  printf("\t\t                      Results table for Re_Profit\n\n");
  printf("\t\t   k     value            e             p             z\n\n");    
  seq( printf( "\t\t%4d%14.8f%14.8f%14.8f%14.8f\n",
               kVals[j],
               Re_Profit[j][1],
               rhs~(Re_Profit[j][2])[]),
       j=1..numelems(kVals)
     );
  printf("\n\n\t\t                      Results table for Int_Profit\n\n");
  printf("\t\t   k     value            e             p             z\n\n");    
  seq( printf( "\t\t%4d%14.8f%14.8f%14.8f%14.8f\n",
               kVals[j],
               Int_Profit[j][1],
               rhs~(Int_Profit[j][2])[]
             ),
       j=1..numelems(kVals)
     );
  printf("\n\n\t\t                       Results table for Range of t\n\n");
  printf( "\t\t                       k=1"):
  printf( "                   k=2"):
  printf( "                   k=3\n\n"):
  printf( "\t\t   t-Range    %9.6f..%9.6f  %9.6f..%9.6f  %9.6f..%9.6f\n",
          seq( [qd_opt[j], t[j]][],
               j=1..numelems(kVals)
             )
        ):
  printf("\n\n\t\t   Results table for Range of t\n\n");
  printf("\t\t   k          range\n\n");
    seq( printf( "\t\t%4d  %9.6f..%9.6f\n",
               kVals[j],
               qd_opt[j],
               t[j]
             ),
       j=1..numelems(kVals)
     );



                                 k=1           k=2           k=3


                      qd_opt   14.10261511   15.36561004   18.03125117
                   Re_Profit  118.99453246  129.08114340  150.27134063
                  Int_Profit  212.61603585  230.51256895  268.05834111
                Whole_Profit  189.50760803  205.90919359  240.42759650
                           t   18.72430066   20.28628510   23.55740010


                                      Results table for Re_Profit

                   k     value            e             p             z

                   1  118.99453246    0.78757903   17.87579032   -0.93338328
                   2  129.08114340    1.70759187   18.53795944   -0.97365482
                   3  150.27134063    2.97811819   19.92706137   -1.04898067


                                      Results table for Int_Profit

                   k     value            e             p             z

                   1  212.61603585    1.04932896   15.49328965   -1.44802041
                   2  230.51256895    2.27401511   16.37007557   -1.47950809
                   3  268.05834111    3.96251316   18.20837721   -1.53511015


                                       Results table for Range of t

                                       k=1                   k=2                   k=3

                   t-Range    14.102615..18.724301  15.365610..20.286285  18.031251..23.557400


                   Results table for Range of t

                   k          range

                   1  14.102615..18.724301
                   2  15.365610..20.286285
                   3  18.031251..23.557400

  
> 

 

``

Download optProb4.mw

 

You mean like this???...

tomleslie 13876
December 22 2018
0 0

@cinderella 

And if it isn't what you want, then you had better think very carefully before your next question, because I'm getting really bored writing trivial formatted print statements

Have you tried reading the help page for printf() ???

> 

#
# Initialize and set up necessary functions
#
  restart;
  with(Optimization):
  y := (p, e) -> alpha-beta*p+k*e:
  g := e-> (1/2)*mu*e^2:
  Lambda:= z-> int((z-u)*phi(u), u=0..z):
  PI := unapply( (p-c)*(y(p, e)-z)-(p-v)*Lambda(z)-g(e), [p,e,z]):
  Pi_rd:= unapply( (p-w)*(y(p, e)-z)-(p-v)*Lambda(z)-g(e), [p, e, z]):
> 

#
# Set up k-values and optimizer conditions
#
  kVals:=[1, 2, 3]:
  optConds:=  p = 0 .. 60,
              e = -1000 .. 30,
              z = -10000 .. 25,
              initialpoint = {e = 0, p = 10, z = 0},
              maximize:

  for ind from 1 by 1 to numelems(kVals) do
    #
    # Set up parameters
    #
      pars:= [ k = kVals[ind],
               beta = 2,
               alpha = 50,
               mu = 10,
               c = 5,
               v = 1,
               phi(u) = 1/2,
               w=10
             ]:
      Int_Profit[ind]:= NLPSolve
                        ( eval
                          ( PI(p, e, z),
                            pars
                          ),
                          optConds
                       );
      qopt[ind]:= eval(z+y(p, e), [ Int_Profit[ind][2][], pars[]]);
      Re_Profit[ind]:= NLPSolve
                       ( eval
                         ( Pi_rd(p, e, z),
                           pars
                         ),
                         optConds
                       );
      qd_opt[ind]:= eval(z+y(p, e), [ Re_Profit[ind][2][], pars[]]);
      M_Profit[ind]:= eval((w-c)*qd_opt[ind], pars);
      Whole_Profit[ind]:= M_Profit[ind]+Re_Profit[ind][1];
      t[ind]:= eval((Int_Profit[ind][1]-Re_Profit[ind][1])/(w-c), pars);
  od:
> 

#
# How to access results, eg for Re_Profit.
#
# Others ie Int_Profit
#           Whole_Profit
#           qopt
#           qd_opt
#           t
#
# are accessed similarly
#
  Re_Profit[1];
  Re_Profit[2];
  Re_Profit[3];

[118.994532457413399, [e = HFloat(0.78757902976968), p = HFloat(17.875790315337227), z = HFloat(-0.9333832842273811)]]

  

[129.081143396166169, [e = HFloat(1.7075918734530875), p = HFloat(18.537959441958364), z = HFloat(-0.9736548233351248)]]

  

[150.271340626426877, [e = HFloat(2.9781181923084055), p = HFloat(19.9270613673776), z = HFloat(-1.0489806678261728)]]

 (1)
> 

  printf( "\n\n\t\t                 k=1           k=2           k=3\n\n");
  printf( "\t\t      qd_opt%14.8f%14.8f%14.8f\n",
          seq( qd_opt[j],
               j=1..numelems(kVals)
             )
        ):
  printf( "\t\t   Re_Profit%14.8f%14.8f%14.8f\n",
          seq( Re_Profit[j][1],
               j=1..numelems(kVals)
             )
        ):
  printf( "\t\t  Int_Profit%14.8f%14.8f%14.8f\n",
          seq( Int_Profit[j][1],
               j=1..numelems(kVals)
             )
        ):
  printf( "\t\tWhole_Profit%14.8f%14.8f%14.8f\n",
          seq( Whole_Profit[j],
               j=1..numelems(kVals)
             )
        ):
  printf( "\t\t           t%14.8f%14.8f%14.8f\n\n\n",
          seq( t[j],
               j=1..numelems(kVals)
             )
        );
  printf("\t\t                      Results table for Re_Profit\n\n");
  printf("\t\t   k     value            e             p             z\n\n");    
  seq( printf( "\t\t%4d%14.8f%14.8f%14.8f%14.8f\n",
               kVals[j],
               Re_Profit[j][1],
               rhs~(Re_Profit[j][2])[]),
       j=1..numelems(kVals)
     );
  printf("\n\n\t\t                      Results table for Int_Profit\n\n");
  printf("\t\t   k     value            e             p             z\n\n");    
  seq( printf( "\t\t%4d%14.8f%14.8f%14.8f%14.8f\n",
               kVals[j],
               Int_Profit[j][1],
               rhs~(Int_Profit[j][2])[]),
       j=1..numelems(kVals)
     );



                                 k=1           k=2           k=3

                      qd_opt   14.10261511   15.36561004   18.03125117
                   Re_Profit  118.99453246  129.08114340  150.27134063
                  Int_Profit  212.61603585  230.51256895  268.05834111
                Whole_Profit  189.50760803  205.90919359  240.42759650
                           t   18.72430066   20.28628510   23.55740010


                                      Results table for Re_Profit

                   k     value            e             p             z

                   1  118.99453246    0.78757903   17.87579032   -0.93338328
                   2  129.08114340    1.70759187   18.53795944   -0.97365482
                   3  150.27134063    2.97811819   19.92706137   -1.04898067


                                      Results table for Int_Profit

                   k     value            e             p             z

                   1  212.61603585    1.04932896   15.49328965   -1.44802041
                   2  230.51256895    2.27401511   16.37007557   -1.47950809
                   3  268.05834111    3.96251316   18.20837721   -1.53511015

  
> 

 

NULL

Download optProb3.mw

 

 

Use simple printf() statements...

tomleslie 13876
December 22 2018
0 0

@cinderella 

As in the attached where I have appended an execution group

> 

#
# Initialize and set up necessary functions
#
  restart;
  with(Optimization):
  y := (p, e) -> alpha-beta*p+k*e:
  g := e-> (1/2)*mu*e^2:
  Lambda:= z-> int((z-u)*phi(u), u=0..z):
  PI := unapply( (p-c)*(y(p, e)-z)-(p-v)*Lambda(z)-g(e), [p,e,z]):
  Pi_rd:= unapply( (p-w)*(y(p, e)-z)-(p-v)*Lambda(z)-g(e), [p, e, z]):
> 

#
# Set up k-values and optimizer conditions
#
  kVals:=[1, 2, 3]:
  optConds:=  p = 0 .. 60,
              e = -1000 .. 30,
              z = -10000 .. 25,
              initialpoint = {e = 0, p = 10, z = 0},
              maximize:

  for ind from 1 by 1 to numelems(kVals) do
    #
    # Set up parameters
    #
      pars:= [ k = kVals[ind],
               beta = 2,
               alpha = 50,
               mu = 10,
               c = 5,
               v = 1,
               phi(u) = 1/2,
               w=10
             ]:
      Int_Profit[ind]:= NLPSolve
                        ( eval
                          ( PI(p, e, z),
                            pars
                          ),
                          optConds
                       );
      qopt[ind]:= eval(z+y(p, e), [ Int_Profit[ind][2][], pars[]]);
      Re_Profit[ind]:= NLPSolve
                       ( eval
                         ( Pi_rd(p, e, z),
                           pars
                         ),
                         optConds
                       );
      qd_opt[ind]:= eval(z+y(p, e), [ Re_Profit[ind][2][], pars[]]);
      M_Profit[ind]:= eval((w-c)*qd_opt[ind], pars);
      Whole_Profit[ind]:= M_Profit[ind]+Re_Profit[ind][1];
      t[ind]:= eval((Int_Profit[ind][1]-Re_Profit[ind][1])/(w-c), pars);
  od:
> 

#
# How to access results, eg for Re_Profit.
#
# Others ie Int_Profit
#           Whole_Profit
#           qopt
#           qd_opt
#           t
#
# are accessed similarly
#
  Re_Profit[1];
  Re_Profit[2];
  Re_Profit[3];

[118.994532457413399, [e = HFloat(0.78757902976968), p = HFloat(17.875790315337227), z = HFloat(-0.9333832842273811)]]

  

[129.081143396166169, [e = HFloat(1.7075918734530875), p = HFloat(18.537959441958364), z = HFloat(-0.9736548233351248)]]

  

[150.271340626426877, [e = HFloat(2.9781181923084055), p = HFloat(19.9270613673776), z = HFloat(-1.0489806678261728)]]

 (1)
> 

printf( "\n\n\t\t                 k=1           k=2           k=3\n\n");
printf("\t\t      qd_opt%14.8f%14.8f%14.8f\n", qd_opt[1], qd_opt[2], qd_opt[3]);
printf("\t\t   Re_Profit%14.8f%14.8f%14.8f\n", Re_Profit[1][1], Re_Profit[2][1], Re_Profit[3][1]);
printf("\t\t  Int_Profit%14.8f%14.8f%14.8f\n", Int_Profit[1][1], Int_Profit[2][1], Int_Profit[3][1]);
printf("\t\tWhole_Profit%14.8f%14.8f%14.8f\n", Whole_Profit[1], Whole_Profit[2], Whole_Profit[3]);
printf("\t\t           t%14.8f%14.8f%14.8f\n", t[1], t[2], t[3]);



                                 k=1           k=2           k=3

                      qd_opt   14.10261511   15.36561004   18.03125117
                   Re_Profit  118.99453246  129.08114340  150.27134063
                  Int_Profit  212.61603585  230.51256895  268.05834111
                Whole_Profit  189.50760803  205.90919359  240.42759650
                           t   18.72430066   20.28628510   23.55740010

  

 

Download optProb2.mw

Further to my earlier post...

tomleslie 13876
December 20 2018
0 0

@nm 

Other than Win7/Win10, I really can't see a difference.

version();
kernelopts(version);
interface(version);
Physics:-Version();

 

User Interface: 1362973
         Kernel: 1362973
        Library: 1362973
                            1362973
  Maple 2018.2, X86 64 WINDOWS, Nov 16 2018, Build ID 1362973
Standard Worksheet Interface, Maple 2018.2, Windows 7, November

   16 2018 Build ID 1362973
"C:\Users\TomLeslie\maple\toolbox\2018\Physics Updates\lib\Physi\

  cs Updates.maple", 2018, December 20, 22:34 hours, version in

   the MapleCloud: 263, version installed in this computer: 263

 

 

 

 

No help but...

tomleslie 13876
December 20 2018
0 0

I just installed this on my Maple 2018.2 version.

Took less than 10secs.

I am running 64-bit Win7, but the fact that the update installs so fast/easily at least suggest that the update (itself) is not the issue?

There are four possibilities...

tomleslie 13876
December 19 2018
0 0

@bliengme 

linear x-axis, linear y-axis

linear x-axis, logarithmic y-axis

logarithmic x-axis, linear y-axis

logarithmic x-axis, logarithmic y-axis

Maple handles all four of these as shown in the attached. So what exactly is the plotting requirement which isn't handled by one of them

(NB plots 'render' better in a Mpale worksheet than on this site)

> 

#
# Linear x-axis, linear y-axis
#
  plot( exp(x),
        x=1..10,
        title= "linear x, linear y",
        titlefont=[times, bold, 20]
      );
#
# Linear x-axis, logarithmic y-axis
#
  plots:-logplot( exp(x),
                  x=1..10,
                  title= "linear x, logarithmic y",
                  titlefont=[times, bold, 20]
                );
#
# Logarithmic x-axis, linear y-axis
#
  plots:-semilogplot( exp(x),
                      x=1..10,
                      title= "logarithmic x, linear y",
                      titlefont=[times, bold, 20]
                    );
#
# Logarithmic x-axis, logarithmic y-axis
#
  plots:-loglogplot( exp(x),
                     x=1..10,
                     title= "logarithmic x, logarithmic y",
                     titlefont=[times, bold, 20]
                   );

 

 

 

 
> 

 

Download logplots.mw

Help yourself...

tomleslie 13876
December 19 2018
0 0

by posting actual code using the big green up-arrow in the Mapleprimes toolbar

If you post a 'picture' of code then you are assuming that someone here is going to retype your code in an attempt to duplicate your problem. Such retyping will rarely (if ever) happen

Observation...

tomleslie 13876
December 19 2018
0 0

@bliengme 

the plots:-logplot() command which provides a linear x-axis and a logarithmic y-axis has been available since Maple VR2 which was released in November 1992!

Just how old is your Maple installation??!!

Might be a version issue??...

tomleslie 13876
December 19 2018
0 0

@siddikmaple 

The error mesage you are getting suggest a problem somewhere in the processing of results, before sending these to Excel files. Since everthing works for me in Maple 2018, this *might* be due to the version of Maple you are running. What is your version?

Try re-executing the attached, then save the results. When prompted for whether you want to save with data, say "yes". Then post the resulting worksheet here

> 

restart
> 

#
# Load packages "properly" and check Maple
# version
#
  with(DirectSearch):
  with(LinearAlgebra):
  kernelopts(version)

`Maple 2018.2, X86 64 WINDOWS, Nov 16 2018, Build ID 1362973`

 (1)
> 

lambda[0] := 0.1e-5

0.1e-5

 (2)
> 

lambda[1] := 0.1e-5

0.1e-5

 (3)
> 

lambda[2] := 0.1e-5

0.1e-5

 (4)
> 

lambda[3] := 0.1e-5

0.1e-5

 (5)
> 

for n to 20 do e1 := p[s] = (n*t[0]/(1-t[0])+n*t[1]/(1-t[1])+n*t[2]/(1-t[2])+n*t[3]/(1-t[3]))*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n; e2 := p[b] = 1-((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n; e3 := r[0] = lambda[0]*(20*p[0, p]*(1-r[0])*((7/2)*p[0, b]-(7/2)*p[0, b]*p[0, c]^5+(15/2)*p[0, c]+(15/2)*p[0, c]^2+(15/2)*p[0, c]^3+(15/2)*p[0, c]^4-30*p[0, c]^5)+20*r[0]*(7/2+(15/2)*p[0, c]+(15/2)*p[0, c]^2+(15/2)*p[0, c]^3+(15/2)*p[0, c]^4-(67/2)*p[0, c]^5)+p[0, b]*(11258.065*p[s]/p[b]+469.725)*(p[0, p]*(1-r[0])*((7/2)*p[0, b]-(7/2)*p[0, b]*p[0, c]^5+(15/2)*p[0, c]+(15/2)*p[0, c]^2+(15/2)*p[0, c]^3+(15/2)*p[0, c]^4-30*p[0, c]^5)+r[0]*(7/2+(15/2)*p[0, c]+(15/2)*p[0, c]^2+(15/2)*p[0, c]^3+(15/2)*p[0, c]^4-(67/2)*p[0, c]^5))+(469.725*(-p[0, p]*r[0]+p[0, p]+r[0]))*(-4*p[0, c]^5+p[0, c]^4+p[0, c]^3+p[0, c]^2+p[0, c])+11727.79-11727.79*p[0, c]^5*(-p[0, p]*r[0]+p[0, p]+r[0])+p[0, c]^5*(20*p[0, p]*(1-r[0])*(7*p[0, b]*(1/2)+30)+3018.625*r[0]+(20*p[0, p]*(1-r[0])*(7*p[0, b]*(1/2)+30)+670*r[0])*p[0, b]*(11258.065*p[s]+469.725*p[b])/(20*p[b])+2348.625*p[0, p]*(1-r[0]))); e4 := r[1] = lambda[1]*(20*p[1, p]*(1-r[1])*((15/2)*p[1, b]-(15/2)*p[1, b]*p[1, c]^5+(31/2)*p[1, c]+(31/2)*p[1, c]^2+(31/2)*p[1, c]^3+(31/2)*p[1, c]^4-62*p[1, c]^5)+20*r[1]*(15/2+(31/2)*p[1, c]+(31/2)*p[1, c]^2+(31/2)*p[1, c]^3+(31/2)*p[1, c]^4-(139/2)*p[1, c]^5)+p[1, b]*(11258.065*p[s]/p[b]+469.725)*(p[1, p]*(1-r[1])*((15/2)*p[1, b]-(15/2)*p[1, b]*p[1, c]^5+(31/2)*p[1, c]+(31/2)*p[1, c]^2+(31/2)*p[1, c]^3+(31/2)*p[1, c]^4-62*p[1, c]^5)+r[1]*(15/2+(31/2)*p[1, c]+(31/2)*p[1, c]^2+(31/2)*p[1, c]^3+(31/2)*p[1, c]^4-(139/2)*p[1, c]^5))+(469.725*(-p[1, p]*r[1]+p[1, p]+r[1]))*(-4*p[1, c]^5+p[1, c]^4+p[1, c]^3+p[1, c]^2+p[1, c])+11727.79-11727.79*p[1, c]^5*(-p[1, p]*r[1]+p[1, p]+r[1])+p[1, c]^5*(20*p[1, p]*(1-r[1])*(15*p[1, b]*(1/2)+62)+3738.625*r[1]+(20*p[1, p]*(1-r[1])*(15*p[1, b]*(1/2)+62)+1390*r[1])*p[1, b]*(11258.065*p[s]+469.725*p[b])/(20*p[b])+2348.625*p[1, p]*(1-r[1]))); e5 := r[2] = lambda[2]*(20*p[2, p]*(1-r[2])*((31/2)*p[2, b]-(31/2)*p[2, b]*p[2, c]^5+(63/2)*p[2, c]+(127/2)*p[2, c]^2+(255/2)*p[2, c]^3+(511/2)*p[2, c]^4-478*p[2, c]^5)+20*r[2]*(31/2+(63/2)*p[2, c]+(127/2)*p[2, c]^2+(255/2)*p[2, c]^3+(511/2)*p[2, c]^4-(987/2)*p[2, c]^5)+p[2, b]*(11258.065*p[s]/p[b]+469.725)*(p[2, p]*(1-r[2])*((31/2)*p[2, b]-(31/2)*p[2, b]*p[2, c]^5+(63/2)*p[2, c]+(127/2)*p[2, c]^2+(255/2)*p[2, c]^3+(511/2)*p[2, c]^4-478*p[2, c]^5)+r[2]*(31/2+(63/2)*p[2, c]+(127/2)*p[2, c]^2+(255/2)*p[2, c]^3+(511/2)*p[2, c]^4-(987/2)*p[2, c]^5))+(469.725*(-p[2, p]*r[2]+p[2, p]+r[2]))*(-4*p[2, c]^5+p[2, c]^4+p[2, c]^3+p[2, c]^2+p[2, c])+11727.79-11727.79*p[2, c]^5*(-p[2, p]*r[2]+p[2, p]+r[2])+p[2, c]^5*(20*p[2, p]*(1-r[2])*(31*p[2, b]*(1/2)+478)+12218.625*r[2]+(20*p[2, p]*(1-r[2])*(31*p[2, b]*(1/2)+478)+9870*r[2])*p[2, b]*(11258.065*p[s]+469.725*p[b])/(20*p[b])+2348.625*p[2, p]*(1-r[2]))); e6 := r[3] = lambda[3]*(20*p[3, p]*(1-r[3])*((31/2)*p[3, b]-(31/2)*p[3, b]*p[3, c]^5+(63/2)*p[3, c]+(127/2)*p[3, c]^2+(255/2)*p[3, c]^3+(511/2)*p[3, c]^4-478*p[3, c]^5)+20*r[3]*(31/2+(63/2)*p[3, c]+(127/2)*p[3, c]^2+(255/2)*p[3, c]^3+(511/2)*p[3, c]^4-(987/2)*p[3, c]^5)+p[3, b]*(11258.065*p[s]/p[b]+469.725)*(p[3, p]*(1-r[3])*((31/2)*p[3, b]-(31/2)*p[3, b]*p[3, c]^5+(63/2)*p[3, c]+(127/2)*p[3, c]^2+(255/2)*p[3, c]^3+(511/2)*p[3, c]^4-478*p[3, c]^5)+r[3]*(31/2+(63/2)*p[3, c]+(127/2)*p[3, c]^2+(255/2)*p[3, c]^3+(511/2)*p[3, c]^4-(987/2)*p[3, c]^5))+(469.725*(-p[3, p]*r[3]+p[3, p]+r[3]))*(-4*p[3, c]^5+p[3, c]^4+p[3, c]^3+p[3, c]^2+p[3, c])+11727.79-11727.79*p[3, c]^5*(-p[3, p]*r[3]+p[3, p]+r[3])+p[3, c]^5*(20*p[3, p]*(1-r[3])*(31*p[3, b]*(1/2)+478)+12218.625*r[3]+(20*p[3, p]*(1-r[3])*(31*p[3, b]*(1/2)+478)+9870*r[3])*p[3, b]*(11258.065*p[s]+469.725*p[b])/(20*p[b])+2348.625*p[3, p]*(1-r[3]))); e7 := t[0] = (-p[0, c]^5+1)/((1-p[0, c])*(r[0]+(1-r[0])/p[0, p]+(-p[0, c]^5+1)/(1-p[0, c])+(7*(p[0, b]*(1-r[0])+r[0]))/(2-2*p[0, b])+(15*p[0, c]^4+15*p[0, c]^3+15*p[0, c]^2+15*p[0, c])/(2-2*p[0, b]))); e8 := t[1] = (-p[1, c]^5+1)/((1-p[1, c])*(r[1]+(1-r[1])/p[1, p]+(-p[1, c]^5+1)/(1-p[1, c])+(15*(p[1, b]*(1-r[1])+r[1]))/(2-2*p[1, b])+(31*p[1, c]^4+31*p[1, c]^3+31*p[1, c]^2+31*p[1, c])/(2-2*p[1, b]))); e9 := t[2] = (-p[2, c]^5+1)/((1-p[2, c])*(r[2]+(1-r[2])/p[2, p]+(-p[2, c]^5+1)/(1-p[2, c])+(31*(p[2, b]*(1-r[2])+r[2]))/(2-2*p[2, b])+(511*p[2, c]^4+255*p[2, c]^3+127*p[2, c]^2+63*p[2, c])/(2-2*p[2, b]))); e10 := t[3] = (-p[3, c]^5+1)/((1-p[3, c])*(r[3]+(1-r[3])/p[3, p]+(-p[3, c]^5+1)/(1-p[3, c])+(31*(p[3, b]*(1-r[3])+r[3]))/(2-2*p[3, b])+(511*p[3, c]^4+255*p[3, c]^3+127*p[3, c]^2+63*p[3, c])/(2-2*p[3, b]))); e11 := p[0, p] = 1-exp(-lambda[0]*(-449.725*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+(11258.065*(n*t[0]/(1-t[0])+n*t[1]/(1-t[1])+n*t[2]/(1-t[2])+n*t[3]/(1-t[3])))*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+469.725)); e12 := p[1, p] = 1-exp(-lambda[1]*(-449.725*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+(11258.065*(n*t[0]/(1-t[0])+n*t[1]/(1-t[1])+n*t[2]/(1-t[2])+n*t[3]/(1-t[3])))*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+469.725)); e13 := p[2, p] = 1-exp(-lambda[2]*(-449.725*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+(11258.065*(n*t[0]/(1-t[0])+n*t[1]/(1-t[1])+n*t[2]/(1-t[2])+n*t[3]/(1-t[3])))*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+469.725)); e14 := p[3, p] = 1-exp(-lambda[3]*(-449.725*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+(11258.065*(n*t[0]/(1-t[0])+n*t[1]/(1-t[1])+n*t[2]/(1-t[2])+n*t[3]/(1-t[3])))*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+469.725)); e15 := p[0, c] = 1-((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[0]); e16 := p[1, c] = 1-((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[1]); e17 := p[2, c] = 1-((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[2]); e18 := p[3, c] = 1-((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[3]); e19 := p[0, b] = 1-(((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[0]))^1; e20 := p[1, b] = 1-(((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[1]))^1; e21 := p[2, b] = 1-(((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[2]))^2; e22 := p[3, b] = 1-(((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[3]))^6; e23 := p[0, s] = n*t[0]*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[0]); e24 := p[1, s] = n*t[1]*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[1]); e25 := p[2, s] = n*t[2]*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[2]); e26 := p[3, s] = n*t[3]*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[3]); CSTR := `~`[`=`]({p[b], p[s], p[0, b], p[0, c], p[0, p], p[0, s], p[1, b], p[1, c], p[1, p], p[1, s], p[2, b], p[2, c], p[2, p], p[2, s], p[3, b], p[3, c], p[3, p], p[3, s], r[0], r[1], r[2], r[3], t[0], t[1], t[2], t[3]}, 0 .. 1); residuals := `~`[lhs-rhs]([e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, e12, e13, e14, e15, e16, e17, e18, e19, e20, e21, e22, e23, e24, e25, e26]); ansLS[n] := Optimization:-LSSolve(residuals, op(CSTR), iterationlimit = 10000); ansDS[n] := DirectSearch:-SolveEquations(residuals, CSTR, initialpoint = ansLS[n][2]) end do
> 

#
# Quick, dirty check on what is being generated by
# above loop - should output two 20x26 matrices
#
  Matrix(20,26,[seq( ansDS[j][3], j=1..10)]);
  Matrix(20,26,[seq( ansLS[j][2], j=1..20)]);

RTABLE(18446744074192895998, anything, Matrix, rectangular, Fortran_order, [], 2, 1 .. 20, 1 .. 26)

  

RTABLE(18446744074192896150, anything, Matrix, rectangular, Fortran_order, [], 2, 1 .. 20, 1 .. 26)

 (6)
> 

#
# Extract the parameter values from the DirectSearch and
# LSSolve results tables. Put these into matrices, together
# with an "index" row and an "index" column.
#
  LSMat:= Matrix( max([indices(ansLS, 'nolist')])+1,
                  numelems(ansLS[1][2])+1
                ):
  LSMat[1,1]:= "n":
  LSMat[1,2..-1]:= convert~( < lhs~(ansLS[1][2])[] >,
                             string
                           ):
  LSMat[2..-1,1]:= < seq
                     ( j,
                       j=1..n-1
                     )
                   >:
  LSMat[2..-1, 2..-1]:= < seq
                          ( Transpose
                            ( < rhs~(ansLS[j][2]) > ),
                            j=1..n-1
                          )
                        >:
  DSMat:= Matrix
          ( max
            ( [ indices
                ( ansDS, 'nolist')
              ]
            )+1,
            numelems
            ( ansDS[1][3] )+1
          ):
  DSMat[1,1]:= "n":
  DSMat[1,2..-1]:= convert~( < lhs~(ansDS[1][3])[] >,
                             string
                           ):
  DSMat[2..-1,1]:= < seq
                     ( j,
                       j=1..n-1
                     )
                   >:
  DSMat[2..-1, 2..-1]:= < seq
                          ( Transpose
                            ( < rhs~(ansDS[j][3]) > ),
                            j=1..n-1
                          )
                        >:
 
> 

#
# Quick, dirty check on what is being generated by
# this code section -should output two 21x27 matrices
#
  LSMat;
  DSMat;

RTABLE(18446744074192896302, anything, Matrix, rectangular, Fortran_order, [], 2, 1 .. 21, 1 .. 27)

  

RTABLE(18446744074192896454, anything, Matrix, rectangular, Fortran_order, [], 2, 1 .. 21, 1 .. 27)

 (7)
> 

#
# So are the results from the LinearSolve() and
# DirectSearch() processes very different?? They
# shouldn't be! For amusement, find the biggest
# absolute difference for any parameter
#
   mi:= max[index](DSMat-LSMat);
  (DSMat-LSMat)[mi];
#
# Differences probably too small to be noticeable
# when data is "rounded" for export to Excel
#

21, 6

  

HFloat(2.2644310899815057e-7)

 (8)
> 

# Send data to two different Excel spreadsheets,
# one for the LSSolve() results and one for the
# DirectSearch() results.
#
# OP will have to change 'fname' to something
# appropriate for his/her machine, OS
#
  fname:= "C:/Users/TomLeslie/Desktop/":
  #fname:= "F:/";
  ExcelTools:-Export( LSMat, cat(fname, "OPTLSres.xlsx"), 1, "B2"):
  ExcelTools:-Export( DSMat, cat(fname, "OPTDSres.xlsx"), 1, "B2"):
> 

 


 

Download toXL2.mw

First 77 78 79 80 81 82 83 Last Page 79 of 207