tomleslie

5194 Reputation

15 Badges

9 years, 244 days

MaplePrimes Activity


These are replies submitted by tomleslie

@radaar 

There are two basic ways to "speed-up" DirectSearch() execution times

  1. For any given set/list of values of the independent variables, reduce the time taken to calculate the objective function. Assume that this objective function will be calculated ~10000 times (at least). Hence saving 0.1sec per function evaluation will save you about 1000secs (or 15 minutes). So spend a lot of time making sure that the objective function evaluation is as fast as possible - this may take quite a lot of thought and careful coding: tough - deal with it
  2. Ensure that the search space of the independent variables is as small as possible. With no constraints on an independent variable, any value between -infinity and infinity is equally possible. Knowing that a value has to be positive will restrict this search space to 0..infinity and you have just reduced the search space by a factor of two (and the execution time by soomething comparable). Knowing that a particular valiable has to be in the range 0..1 reduces the search space dramatically, with hugely beneficial effects on execution time

For "abstract mathematical" optimisation problems search space reduction may be inappropriate. For any "real-world" problems, the "real world" will set constraints on possible values for the variables - use this information ad set the constraints as tightly as possible.

@David Sycamore 

The attached code runs in

Maple 18
Maple 2015
Maple 2016
Maple 2017
Maple 2018
Maple 2019

I know this because I have just checked it in all of these releases. The version uploaded here was run in  Maple 2017 which you claim to have, This can be verified from the output of the interface(version) command, in the following, which is

`Standard Worksheet Interface, Maple 2017.3, Windows 7, September 13 2017 Build ID 1262472`

(1)

I have added a couple of execution groups at the end, just to run the main working procedure on the numbers obtained by Carl Love (more out of interest, because I'm pretty sure he is correct!).

If you find that this code *still* does not execute correctly, then I suggest you perform the following

  1. download the code using the link at the bottom of this post
  2. execute the code - just use the !!! toolbar button in Maple
  3. save the resulting worksheet (any filename will do)
  4. upload this saved worksheet here using the big green up-arrow in the Mapleprimes toolbar

  restart;
  interface(version);
#
# Define the upper limit for the odd numbers
# to be tested and the upper limit for the
# even addends
#
  Nodd:=10000:
  Neven:=10000:
#
# Procedure which does the work
#
  doWork:= proc( p::posint, N::posint )
                 local #
                       # Edited this to run in Maple 18
                       # whihc means it will probably run
                       # in every subsequent version
                       #
                       # f:=NumberTheory:-PrimeFactors(p),
                        f:=numtheory:-factorset(p),
                       j:
                 for j from 2 by 2 to N do
                     if   isprime~(f+~j)={true}
                     then #
                          # Return a list containing the
                          # input odd number, its prime
                          # factors, the lowest even number
                          # which produces a a new set of primes,
                          # and this new set of primes
                          #
                            return [p, f, j, f+~j];
                     fi;
                 od;
                 return #
                        # No even addend found: return a list
                        # containing the input odd number and
                        # a zero
                        #
                          [p, f, 0];
           end proc:

`Standard Worksheet Interface, Maple 2017.3, Windows 7, September 13 2017 Build ID 1262472`

(1)

#
# Check how long it takes to return the complete list
# of answers for the specific values of NO and NE
#
# Results on my machine were. NB all timings approximate
#
#    Nodd    Neven      real time
#
#    100     100        15  msecs
#    100     1000       15  msecs
#    100     10000      15  msecs
#   
#    1000    100        66  msecs
#    1000    1000       77  msecs
#    1000    10000      1.0 secs
#  
#    10000   100        0.5 secs
#    10000   1000       1.5 secs
#    10000   10000      16  secs     
#
  ans:=CodeTools:-Usage([seq( doWork(j, Neven), j=3..Nodd,2)]):

memory used=1.73GiB, alloc change=40.00MiB, cpu time=14.27s, real time=14.27s, gc time=655.20ms

 

#
# Check answers for a couple of the cases used by OP
# to illustrate problem
#
# Utility to facilitate lookup of the answer in the
# main result list for any supplied odd number
#
  getVal:=p->(p-1)/2:
  ans[ getVal(119) ];
  ans[ getVal(105) ];

[119, {7, 17}, 6, {13, 23}]

 

[105, {3, 5, 7}, 0]

(2)

#
# Split the original list into successes and failures.
# Out of idle curiosity, check the number of successes
# and failures
#
  success,failure:=selectremove(i->numelems(i)=4, ans):
  numelems(success);
  numelems(failure);
#
# Output the first few successes and failures
#
  success[1..100];
  failure[1..100];
#
# Check whether the failure cases *always* have '3' as
# a prime factor. This should output any odd number
# which "failed" to generate a new set of prime numbers,
# and itself doesn't have '3' as a prime factor.
#
# It outputs nothing  interesting? coincidence??
#
  seq( `if`( j[2][1]=3, NULL, j[1]), j in failure);

4534

 

465

 

[[3, {3}, 2, {5}], [5, {5}, 2, {7}], [7, {7}, 4, {11}], [9, {3}, 2, {5}], [11, {11}, 2, {13}], [13, {13}, 4, {17}], [15, {3, 5}, 2, {5, 7}], [17, {17}, 2, {19}], [19, {19}, 4, {23}], [21, {3, 7}, 4, {7, 11}], [23, {23}, 6, {29}], [25, {5}, 2, {7}], [27, {3}, 2, {5}], [29, {29}, 2, {31}], [31, {31}, 6, {37}], [33, {3, 11}, 2, {5, 13}], [35, {5, 7}, 6, {11, 13}], [37, {37}, 4, {41}], [39, {3, 13}, 4, {7, 17}], [41, {41}, 2, {43}], [43, {43}, 4, {47}], [45, {3, 5}, 2, {5, 7}], [47, {47}, 6, {53}], [49, {7}, 4, {11}], [51, {3, 17}, 2, {5, 19}], [53, {53}, 6, {59}], [55, {5, 11}, 2, {7, 13}], [57, {3, 19}, 4, {7, 23}], [59, {59}, 2, {61}], [61, {61}, 6, {67}], [63, {3, 7}, 4, {7, 11}], [65, {5, 13}, 6, {11, 19}], [67, {67}, 4, {71}], [69, {3, 23}, 8, {11, 31}], [71, {71}, 2, {73}], [73, {73}, 6, {79}], [75, {3, 5}, 2, {5, 7}], [77, {7, 11}, 6, {13, 17}], [79, {79}, 4, {83}], [81, {3}, 2, {5}], [83, {83}, 6, {89}], [85, {5, 17}, 2, {7, 19}], [87, {3, 29}, 2, {5, 31}], [89, {89}, 8, {97}], [91, {7, 13}, 4, {11, 17}], [93, {3, 31}, 10, {13, 41}], [95, {5, 19}, 12, {17, 31}], [97, {97}, 4, {101}], [99, {3, 11}, 2, {5, 13}], [101, {101}, 2, {103}], [103, {103}, 4, {107}], [107, {107}, 2, {109}], [109, {109}, 4, {113}], [111, {3, 37}, 4, {7, 41}], [113, {113}, 14, {127}], [115, {5, 23}, 6, {11, 29}], [117, {3, 13}, 4, {7, 17}], [119, {7, 17}, 6, {13, 23}], [121, {11}, 2, {13}], [123, {3, 41}, 2, {5, 43}], [125, {5}, 2, {7}], [127, {127}, 4, {131}], [129, {3, 43}, 4, {7, 47}], [131, {131}, 6, {137}], [133, {7, 19}, 4, {11, 23}], [135, {3, 5}, 2, {5, 7}], [137, {137}, 2, {139}], [139, {139}, 10, {149}], [141, {3, 47}, 14, {17, 61}], [143, {11, 13}, 6, {17, 19}], [145, {5, 29}, 2, {7, 31}], [147, {3, 7}, 4, {7, 11}], [149, {149}, 2, {151}], [151, {151}, 6, {157}], [153, {3, 17}, 2, {5, 19}], [155, {5, 31}, 6, {11, 37}], [157, {157}, 6, {163}], [159, {3, 53}, 8, {11, 61}], [161, {7, 23}, 6, {13, 29}], [163, {163}, 4, {167}], [165, {3, 5, 11}, 2, {5, 7, 13}], [167, {167}, 6, {173}], [169, {13}, 4, {17}], [171, {3, 19}, 4, {7, 23}], [173, {173}, 6, {179}], [175, {5, 7}, 6, {11, 13}], [177, {3, 59}, 2, {5, 61}], [179, {179}, 2, {181}], [181, {181}, 10, {191}], [183, {3, 61}, 10, {13, 71}], [185, {5, 37}, 6, {11, 43}], [187, {11, 17}, 2, {13, 19}], [189, {3, 7}, 4, {7, 11}], [191, {191}, 2, {193}], [193, {193}, 4, {197}], [197, {197}, 2, {199}], [199, {199}, 12, {211}], [201, {3, 67}, 4, {7, 71}], [203, {7, 29}, 12, {19, 41}], [205, {5, 41}, 2, {7, 43}]]

 

[[105, {3, 5, 7}, 0], [195, {3, 5, 13}, 0], [231, {3, 7, 11}, 0], [285, {3, 5, 19}, 0], [315, {3, 5, 7}, 0], [357, {3, 7, 17}, 0], [429, {3, 11, 13}, 0], [465, {3, 5, 31}, 0], [483, {3, 7, 23}, 0], [525, {3, 5, 7}, 0], [555, {3, 5, 37}, 0], [585, {3, 5, 13}, 0], [609, {3, 7, 29}, 0], [627, {3, 11, 19}, 0], [645, {3, 5, 43}, 0], [663, {3, 13, 17}, 0], [693, {3, 7, 11}, 0], [735, {3, 5, 7}, 0], [855, {3, 5, 19}, 0], [861, {3, 7, 41}, 0], [897, {3, 13, 23}, 0], [915, {3, 5, 61}, 0], [945, {3, 5, 7}, 0], [969, {3, 17, 19}, 0], [975, {3, 5, 13}, 0], [987, {3, 7, 47}, 0], [1005, {3, 5, 67}, 0], [1023, {3, 11, 31}, 0], [1071, {3, 7, 17}, 0], [1095, {3, 5, 73}, 0], [1113, {3, 7, 53}, 0], [1131, {3, 13, 29}, 0], [1155, {3, 5, 7, 11}, 0], [1185, {3, 5, 79}, 0], [1221, {3, 11, 37}, 0], [1239, {3, 7, 59}, 0], [1287, {3, 11, 13}, 0], [1311, {3, 19, 23}, 0], [1365, {3, 5, 7, 13}, 0], [1395, {3, 5, 31}, 0], [1419, {3, 11, 43}, 0], [1425, {3, 5, 19}, 0], [1449, {3, 7, 23}, 0], [1455, {3, 5, 97}, 0], [1491, {3, 7, 71}, 0], [1545, {3, 5, 103}, 0], [1575, {3, 5, 7}, 0], [1581, {3, 17, 31}, 0], [1599, {3, 13, 41}, 0], [1617, {3, 7, 11}, 0], [1635, {3, 5, 109}, 0], [1653, {3, 19, 29}, 0], [1665, {3, 5, 37}, 0], [1743, {3, 7, 83}, 0], [1755, {3, 5, 13}, 0], [1785, {3, 5, 7, 17}, 0], [1827, {3, 7, 29}, 0], [1833, {3, 13, 47}, 0], [1869, {3, 7, 89}, 0], [1881, {3, 11, 19}, 0], [1887, {3, 17, 37}, 0], [1905, {3, 5, 127}, 0], [1935, {3, 5, 43}, 0], [1989, {3, 13, 17}, 0], [1995, {3, 5, 7, 19}, 0], [2013, {3, 11, 61}, 0], [2067, {3, 13, 53}, 0], [2079, {3, 7, 11}, 0], [2085, {3, 5, 139}, 0], [2121, {3, 7, 101}, 0], [2139, {3, 23, 31}, 0], [2145, {3, 5, 11, 13}, 0], [2193, {3, 17, 43}, 0], [2205, {3, 5, 7}, 0], [2211, {3, 11, 67}, 0], [2247, {3, 7, 107}, 0], [2265, {3, 5, 151}, 0], [2301, {3, 13, 59}, 0], [2325, {3, 5, 31}, 0], [2337, {3, 19, 41}, 0], [2355, {3, 5, 157}, 0], [2373, {3, 7, 113}, 0], [2409, {3, 11, 73}, 0], [2415, {3, 5, 7, 23}, 0], [2445, {3, 5, 163}, 0], [2499, {3, 7, 17}, 0], [2535, {3, 5, 13}, 0], [2541, {3, 7, 11}, 0], [2553, {3, 23, 37}, 0], [2565, {3, 5, 19}, 0], [2583, {3, 7, 41}, 0], [2607, {3, 11, 79}, 0], [2625, {3, 5, 7}, 0], [2679, {3, 19, 47}, 0], [2691, {3, 13, 23}, 0], [2697, {3, 29, 31}, 0], [2715, {3, 5, 181}, 0], [2745, {3, 5, 61}, 0], [2751, {3, 7, 131}, 0], [2769, {3, 13, 71}, 0]]

(3)

#
# What happens when the procedure doWork is applied to the
# odd numbers found by Carl Love. Looks like he was right.
# Final value in each output list is 0, indicating that no even
# number (up to 10000) added to the factor list produces a
# a list of primes.
#
# I suppose there *might* be even numbers greater
# than the tested limit (10000) which could prove the
# opposite
#
  doWork~( [ 105, 95095, 215656441, 5070519693819151, 121621899527020215557,
             21664759328930363899339259137, 6118163487924060025029213220329793,
             3042203271558217610793849843552106188681707,
             410050702193824620789441168165116752238670790985538614407,
             243521388535217606004348248409375049907078159488916219637351331,
             65167879497494548612895080562902666093018940639437608943583983945258079733989
           ],
           10000
         );

[[105, {3, 5, 7}, 0], [95095, {5, 7, 11, 13, 19}, 0], [215656441, {7, 11, 13, 17, 19, 23, 29}, 0], [5070519693819151, {11, 13, 17, 19, 23, 29, 31, 37, 43, 47, 71}, 0], [121621899527020215557, {13, 17, 19, 23, 29, 31, 37, 41, 47, 53, 59, 61, 103}, 0], [21664759328930363899339259137, {17, 19, 23, 29, 31, 37, 41, 43, 47, 59, 61, 67, 73, 79, 83, 89, 103}, 0], [6118163487924060025029213220329793, {19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 71, 73, 83, 89, 101, 103, 131, 163, 191}, 0], [3042203271558217610793849843552106188681707, {23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 97, 101, 103, 109, 127, 131, 137, 157, 233}, 0], [410050702193824620789441168165116752238670790985538614407, {29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 97, 103, 107, 109, 113, 127, 139, 149, 151, 173, 179, 181, 191, 193, 317}, 0], [243521388535217606004348248409375049907078159488916219637351331, {31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 107, 113, 127, 131, 137, 139, 149, 157, 173, 179, 181, 271, 277, 311, 577}, 0], [65167879497494548612895080562902666093018940639437608943583983945258079733989, {37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 131, 137, 139, 149, 151, 167, 173, 179, 193, 197, 199, 229, 233, 239, 277, 317, 383, 439}, 0]]

(4)

doWork(81345, 10000);

[81345, {3, 5, 11, 17, 29}, 2, {5, 7, 13, 19, 31}]

(5)

 

Download prmeProb3.mw

 

 

@David Sycamore 

My original code was developed (and runs) in Maple 2019. However it uses the NumberTheory() package which was introduced in Maple 2016.

Attached is a slightly modified version using the older numtheory() package which runs in Maple 18 (and probably earlier although I can't test because Maple 18 is the oldest version I have installed)

I reccommend that "generall" you don't "copy" the code. Scroll all the way to the bottom of this message and use the blue download link. Your browser should ask you what application you wish to use for this code - Just tell it Maple.

( In the attached I have removed all of the "output" so that you can copy/paste if you absolutely have to.)

You will have to use the code, as far as (and including) the command

ans:=CodeTools:-Usage([seq( doWork(j, Neven), j=3..Nodd,2)]):

which is where the calculation is actually performed, ie where the doWork() procedure is actually executed. The above command runs doWork() for every odd number between 3 and Nodd. The CodeToolsUsage() wrapper on the above command is just for timing purposes.

Commands after this point

  1. make accessing this results list  a little easier
  2. split the single returned list into two - "success" and "failure"
  3. check whether the failure list contains an entry in the which doesn't have three as a prime factor (It doesn't)

  restart;
  interface(version);
#
# Define the upper limit for the odd numbers
# to be tested and the upper limit for the
# even addends
#
  Nodd:=10000:
  Neven:=10000:
#
# Procedure which does the work
#
  doWork:= proc( p::posint, N::posint )
                 local #
                       # Edited this to run in Maple 18
                       #
                       # f:=NumberTheory:-PrimeFactors(p),
                        f:=numtheory:-factorset(p),
                       j:
                 for j from 2 by 2 to N do
                     if   isprime~(f+~j)={true}
                     then #
                          # Return a list containing the
                          # input odd number, its prime
                          # factors, the lowest even number
                          # which produces a a new set of primes,
                          # and this new set of primes
                          #
                            return [p, f, j, f+~j];
                     fi;
                 od;
                 return #
                        # No even addend found: return a list
                        # containing the input odd number and
                        # a zero
                        #
                          [p, f, 0];
           end proc:

#
# Check how long it takes to return the complete list
# of answers for the specific values of NO and NE
#
# Results on my machine were. NB all timings approximate
#
#    Nodd    Neven      real time
#
#    100     100        15  msecs
#    100     1000       15  msecs
#    100     10000      15  msecs
#   
#    1000    100        66  msecs
#    1000    1000       77  msecs
#    1000    10000      1.0 secs
#  
#    10000   100        0.5 secs
#    10000   1000       1.5 secs
#    10000   10000      16  secs     
#
  ans:=CodeTools:-Usage([seq( doWork(j, Neven), j=3..Nodd,2)]):

#
# Check answers for a couple of the cases used by OP
# to illustrate problem
#
# Utility to facilitate lookup of the answer in the
# main result list for any supplied odd number
#
  getVal:=p->(p-1)/2:
  ans[ getVal(119) ];
  ans[ getVal(105) ];

#
# Split the original list into successes and failures.
# Out of idle curiosity, check the number of successes
# and failures
#
  success,failure:=selectremove(i->numelems(i)=4, ans):
  numelems(success);
  numelems(failure);
#
# Output the first few successes and failures
#
  success[1..100];
  failure[1..100];
#
# Check whether the failure cases *always* have '3' as
# a prime factor. This should output any odd number
# which "failed" to generate a new set of prime numbers,
# and itself doesn't have '3' as a prime factor.
#
# It outputs nothing  interesting? coincidence??
#
  seq( `if`( j[2][1]=3, NULL, j[1]), j in failure);

 

Download prmeProb2.mw

 

 

@ogunmiloro 

The routine getPlot() will accept a variable name, a list of variable values, and the name of the dependent function you want to plot. It can be used in a variety of ways as shown in the attached

  restart;
  interface(version);

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

(1)

#
# "Default" parameter values
#
  params:= { M__h=.50, beta__o=0.34e-1, beta__1=0.25e-1,
             mu__r=0.4e-3, sigma=.7902, alpha=.11,
             psi=0.136e-3, xi=0.5e-1, Ggamma=.7,
             M__c=.636, mu__b=0.5e-2
           }:
#
# ODEs and bcs
#
  ODEs:= diff(B(T), T) = M__h-beta__1*psi*B(T)*G(T)-mu__r*B(T),
         diff(C(T), T) = beta__1*psi*B(T)*G(T)-sigma*psi*beta__1*E(T)*DD(T)-(alpha+xi+mu__r)*C(T),
         diff(DD(T), T) = alpha*C(T)-(Ggamma+mu__r)*DD(T),
         diff(E(T), T) = Ggamma*DD(T)+sigma*psi*beta__1*E(T)*G(T)-mu__r*E(T),
         diff(F(T), T) = M__c-psi*beta__o*F(T)*C(T)-mu__b*F(T),
         diff(G(T), T) = psi*beta__o*F(T)*C(T)-mu__b*G(T):
  bcs:= B(0) = .50, C(0) = .30, DD(0) = .21, E(0) = .14,
        F(0) = .70, G(0) = .45:
#
# Routine which accepts a parameter name, a list of
# parameter values, and the name of a function to plot
#  
  getPlot:= proc( pn::name, pv::list, fv::function )
                 uses plots;
                 local plts:=NULL,
                       colors:=[ red, green, blue, orange,
                                 brown, black, cyan
                               ],
                       k,
                       dSol:= pname-> dsolve
                                      ( [ eval
                                          ( ODEs,
                                            [ seq
                                              ( `if`( lhs(j)=pname, NULL, j),
                                                j in params
                                              )
                                            ]
                                          ),
                                          bcs
                                        ],
                                        numeric,
                                        parameters=[ pname]
                                      ),
                       ans:=dSol(pn):
                 for k from 1 to numelems(pv) do
                     ans(parameters=[pv[k]]):
                     plts:=plts, odeplot
                                 ( ans,
                                   [T, fv],
                                   T = 0 .. 30,
                                   color=colors[k]
                                 );
                 end do:
                 return display( [plts],
                                 title = cat( convert(fv, string),
                                              " as ",
                                              convert(pn, string),
                                              " varies"
                                            ),
                                 titlefont = [times, bold, 20]
                               ):
           end proc:

#
# Usage examples. The arguments are
#
# 1. the name of the parameter you want to vary
# 2. the values for this parameter (don't use
#    more than seven values, becuase only seven
#    colors are defined in the routine getPlot()
# 3. the function you want to be plotted
#
  getPlot( xi, [0.25e-1, 0.5e-1, 0.75e-1], E(T) );
  getPlot( xi, [0.25e-1, 0.5e-1, 0.75e-1], C(T) );
  getPlot( Ggamma, [0.3, 0.7, 1.0], E(T));

 

 

 

#
# Or even (assuming you want plots of all the
# dependent functions)
#
  for j in indets([ODEs], function(name)) do
      getPlot( M__h, [0.3, 0.4, 0.5, 0.6, 0.7], j );
  od;

 

 

 

 

 

 

 


 

Download odeRes4.mw

 

 

@ogunmiloro 

In the attached I have addressed the following problems

  1. Code now runs correctly in Maple 18. My previous worksheet was developed in Maple 2019, where the interface(rtablesize) command has more flexibility. I have redefined this command in a way that ensures the output matrix will appear in Maple 18
  2. I have added an execution group which shows how to genearate plots of all independent variables as one of the system parameters varies. I chose the parameter 'xi' and produced plots for 'xi=[0.25, 0.5,0.75]. No particular reason for selecting this parameter, or the list of values. I just needed something which would illustrate the principle
  3. As I have already stated the "power series" solution whihc you include in your worksheet is pretty much nonsense. I don't care what you think it represents. You know it gives completely incorrect answers. I know it gives completely incorrect answers. You have two choices
    1. The "differential transform method" is complete rubbich (doubtful)
    2. You can't code the differetnial transform method

Pick one of the above two options, becuase you have no other choices. In the attache I have retained the code which shows what the power series solutions for your problem actually are, and why these are nowhere near as good as the default Runge-Kutta method used by the dsolve(..numeric...) command. I still have no idea why you are interested in generating inaccurate power series solutions, when accurated Runge-Kutta solutions are available


 

restart; with(plots); interface(version); _local(gamma)

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

(1)

M__h := .50; beta__o := 0.34e-1; beta__1 := 0.25e-1; mu__r := 0.4e-3; sigma := .7902; alpha := .11; psi := 0.136e-3; xi := 0.5e-1; gamma := .7; M__c := .636; mu__b := 0.5e-2

ODEs := diff(B(T), T) = M__h-beta__1*psi*B(T)*G(T)-mu__r*B(T), diff(C(T), T) = beta__1*psi*B(T)*G(T)-sigma*psi*beta__1*E(T)*DD(T)-(alpha+xi+mu__r)*C(T), diff(DD(T), T) = alpha*C(T)-(gamma+mu__r)*DD(T), diff(E(T), T) = gamma*DD(T)+sigma*psi*beta__1*E(T)*G(T)-mu__r*E(T), diff(F(T), T) = M__c-psi*beta__o*F(T)*C(T)-mu__b*F(T), diff(G(T), T) = psi*beta__o*F(T)*C(T)-mu__b*G(T); bcs := B(0) = .50, C(0) = .30, DD(0) = .21, E(0) = .14, F(0) = .70, G(0) = .45

ans := dsolve([ODEs, bcs], numeric); for j in indets([ODEs], function(name)) do odeplot(ans, [T, j], T = 0 .. 30, title = convert(j, string), titlefont = [tims, bold, 20], thickness = 2, color = red) end do

 

 

 

 

 

 

interface(rtablesize = 22); M := Matrix([`~`[lhs](ans(0)), seq(`~`[rhs](ans(j)), j = 0 .. 20)])

Matrix(%id = 18446744074356196526)

(2)

#
# Generate power series solutions of order 12
# (Maple default is 6) for all ODEs
#
  Order:=12:
  ans2:= convert~
         ( evalf~
           ( dsolve
             ( {ODEs, bcs},
               indets([ODEs], function(name)),
               'series'
             )
           ),
           polynom
         );
#
# Plot the "numeric" solution obtained earlier
# and the equivalent power series solution
# obtaine above on th same graph.
#
# Note that the range of the plots is reduced,
# because even with 12-th order polynomials,
# the power series soluton and the numerical
# solution start to deviate for T>~6
#
  for j in indets([ODEs], function(name)) do
      display
      ( [ odeplot
          ( ans,
            [T, j],
            T = 0 .. 10,
            thickness = 2,
            color = red
          ),
          plot
          ( eval(j, ans2),
            T = 0 .. 10,
            thickness = 2,
            color=blue
          )
        ],
        title = cat(convert(j,string), " (Numerical and PowerSeries compared)"),
        titlefont = [tims, bold, 20]
      )
   end do;

{B(T) = .5000000000+.4997992350*T-0.1003402817e-3*T^2+0.1470042520e-7*T^3-0.4186689729e-11*T^4+0.1797722185e-13*T^5-0.7652371692e-15*T^6+0.2923183302e-16*T^7-0.8844182475e-18*T^8+0.2215403975e-19*T^9-0.4589870880e-21*T^10+0.5389509965e-23*T^11, C(T) = .3000000000-0.4811931399e-1*T+0.3859529419e-2*T^2-0.2063396439e-3*T^3+0.8267055682e-5*T^4-0.2631388485e-6*T^5+0.6554216926e-8*T^6-0.5551714044e-10*T^7-0.1520371564e-10*T^8+0.2777703195e-11*T^9-0.3923000816e-12*T^10+0.4970214040e-13*T^11, DD(T) = .2100000000-.1140840000*T+0.3730565453e-1*T^2-0.8568110732e-2*T^3+0.1494601849e-2*T^4-0.2091819518e-3*T^5+0.2441368229e-4*T^6-0.2442660302e-5*T^7+0.2138541461e-6*T^8-0.1664279070e-7*T^9+0.1165691616e-8*T^10-0.7422668732e-10*T^11, E(T) = .1400000000+.1469441693*T-0.3995870043e-1*T^2+0.8709964152e-2*T^3-0.1500287681e-2*T^4+0.2093639084e-3*T^5-0.2441847491e-4*T^6+0.2442759169e-5*T^7-0.2138545263e-6*T^8+0.1664257448e-7*T^9-0.1165658906e-8*T^10+0.7422262550e-10*T^11, F(T) = .7000000000+.6324990290*T-0.1581608397e-2*T^2+0.2679492283e-5*T^3-0.6093270788e-8*T^4+0.1272029637e-9*T^5-0.4253521584e-11*T^6+0.1189725836e-12*T^7-0.2702521878e-14*T^8+0.3074973933e-16*T^9+0.3481499101e-17*T^10-0.6346661500e-18*T^11, G(T) = .4500000000-0.2249028960e-2*T+0.5983397429e-5*T^2-0.5345061622e-7*T^3+0.2810718704e-8*T^4-0.1239204116e-9*T^5+0.4250786124e-11*T^6-0.1189706297e-12*T^7+0.2702520657e-14*T^8-0.3074973865e-16*T^9-0.3481499102e-17*T^10+0.6346661500e-18*T^11}

 

 

 

 

 

 

 

#
# Start again, and this time set thngs up so that
# one can vary one of the parameters and produce
# corresponding plots of all the independent variables
# as the specified parameer varies
#
  restart;

  with(plots):
  interface(version);

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

(3)

#
# Parameter values as a "reminder". NOte that none
# of these parameter values are assigned
#
# M__h := .50; beta__o := 0.34e-1; beta__1 := 0.25e-1;
# mu__r := 0.4e-3; sigma := .7902; alpha := .11;
# psi := 0.136e-3; xi := 0.5e-1; Ggamma := .7;
# M__c := .636; mu__b := 0.5e-2;

#
# Define the ODEs, ics, and solve the ODE system
# with no determined parameter values
#
  ODEs:= diff(B(T), T) = M__h-beta__1*psi*B(T)*G(T)-mu__r*B(T),
         diff(C(T), T) = beta__1*psi*B(T)*G(T)-sigma*psi*beta__1*E(T)*DD(T)-(alpha+xi+mu__r)*C(T),
         diff(DD(T), T) = alpha*C(T)-(Ggamma+mu__r)*DD(T),
         diff(E(T), T) = Ggamma*DD(T)+sigma*psi*beta__1*E(T)*G(T)-mu__r*E(T),
         diff(F(T), T) = M__c-psi*beta__o*F(T)*C(T)-mu__b*F(T),
         diff(G(T), T) = psi*beta__o*F(T)*C(T)-mu__b*G(T):
  bcs:= B(0) = .50, C(0) = .30, DD(0) = .21, E(0) = .14,
        F(0) = .70, G(0) = .45:
  ans:= dsolve( [ODEs, bcs],
                numeric,
                parameters=[ M__h, beta__o, beta__1,
                             mu__r, sigma, alpha,
                             psi, xi, Ggamma,
                             M__c, mu__b
                           ]
              ):
#
# Define a listlist of parameter values. Note that the
# parameter 'xi' takes the values [0.25, 0.5, 0.75]
# No real reason for picking this particular parameter
# to vary, or the specific values - just need something
# for illustration
#
  pVals:= [ [ .50, 0.34e-1, 0.25e-1, 0.4e-3, .7902, .11,
               0.136e-3, 0.25e-1, .7, .636, 0.5e-2],
            [ .50, 0.34e-1, 0.25e-1, 0.4e-3, .7902, .11,
               0.136e-3, 0.5e-1, .7, .636, 0.5e-2],
            [ .50, 0.34e-1, 0.25e-1, 0.4e-3, .7902, .11,
               0.136e-3, 0.75e-1, .7, .636, 0.5e-2]
          ]:
  colors:=[red, green, blue]:
#
# Generate plots for each dependent variable as the
# parameter 'xi' varies
#
  for j in indets([ODEs], function(name)) do
      plts:=NULL:
      for k from 1 to 3 do
          ans(parameters=pVals[k]):
          plts:= plts,
                 odeplot
                 ( ans,
                   [T, j],
                   T = 0 .. 30,
                   color=colors[k]
                 );
      end do:
      display( [plts],
               title = convert(j, string),
               titlefont = [tims, bold, 20],
               thickness = 2
             );
  end do;

"plts:="

 

 

"plts:="

 

 

"plts:="

 

 

"plts:="

 

 

"plts:="

 

 

"plts:="

 

 

 

 


 

Download odeRes3.mw

 

  1. It is possible to ask Maple to generate "power series" solutions for systems of ODEs. This is done in the attached, where I have generated the relevant series up to Order(12). However as the comparison graphs demonstrate, the power series solution deviates substantially from the numerical solution obtained by Runge-Kutta methods, as T increases. The attached shows that "reasonable" agreement between the two methods is only obtained up to ~6. One could increase the order of the power series solutions (the attached uses 12) and the range of "reasonable" will be extended.
  2. The code you supply in your worksheet does not implement the method required to produce a power series solutiion for an ODE system. In fact I have absolutely no idea what your code calculates, but I guarantee that it has nothing to do with power series solution of the ODE system.

If you really want to develop your own code for producing power series solutions, I suggest you read

  1. Chaper 3 of Elementary Differention Equations by Edwards  and  Penney, or
  2. https://en.wikipedia.org/wiki/Power_series_solution_of_differential_equations (for a quick introduction).

before starting

restart; with(plots); _local(gamma)

M__h := .50; beta__o := 0.34e-1; beta__1 := 0.25e-1; mu__r := 0.4e-3; sigma := .7902; alpha := .11; psi := 0.136e-3; xi := 0.5e-1; gamma := .7; M__c := .636; mu__b := 0.5e-2

ODEs := diff(B(T), T) = M__h-beta__1*psi*B(T)*G(T)-mu__r*B(T), diff(C(T), T) = beta__1*psi*B(T)*G(T)-sigma*psi*beta__1*E(T)*DD(T)-(alpha+xi+mu__r)*C(T), diff(DD(T), T) = alpha*C(T)-(gamma+mu__r)*DD(T), diff(E(T), T) = gamma*DD(T)+sigma*psi*beta__1*E(T)*G(T)-mu__r*E(T), diff(F(T), T) = M__c-psi*beta__o*F(T)*C(T)-mu__b*F(T), diff(G(T), T) = psi*beta__o*F(T)*C(T)-mu__b*G(T); bcs := B(0) = .50, C(0) = .30, DD(0) = .21, E(0) = .14, F(0) = .70, G(0) = .45

ans := dsolve([ODEs, bcs], numeric); for j in indets([ODEs], function(name)) do odeplot(ans, [T, j], T = 0 .. 30, title = convert(j, string), titlefont = [tims, bold, 20], thickness = 2, color = red) end do

 

 

 

 

 

 

interface(rtablesize = [22, 10]); M := Matrix([`~`[lhs](ans(0)), seq(`~`[rhs](ans(j)), j = 0 .. 20)])

Matrix(%id = 18446744074590329062)

(1)

#
# Generate power series solutions of order 12
# (Maple default is 6) for all ODEs
#
  Order:=12:
  ans2:= convert~
         ( evalf~
           ( dsolve
             ( {ODEs, bcs},
               indets([ODEs], function(name)),
               'series'
             )
           ),
           polynom
         );
#
# Plot the "numeric" solution obtained earlier
# and the equivalent power series solution
# obtaine above on th same graph.
#
# Note that the range of the plots is reduced,
# because even with 12-th order polynomials,
# the power series soluton and the numerical
# solution start to deviate for T>~6
#
  for j in indets([ODEs], function(name)) do
      display
      ( [ odeplot
          ( ans,
            [T, j],
            T = 0 .. 10,
            thickness = 2,
            color = red
          ),
          plot
          ( eval(j, ans2),
            T = 0 .. 10,
            thickness = 2,
            color=blue
          )
        ],
        title = cat(convert(j,string), " (Numerical and PowerSeries compared)"),
        titlefont = [tims, bold, 20]
      )
   end do

{B(T) = .5000000000+.4997992350*T-0.1003402817e-3*T^2+0.1470042520e-7*T^3-0.4186689729e-11*T^4+0.1797722185e-13*T^5-0.7652371692e-15*T^6+0.2923183302e-16*T^7-0.8844182475e-18*T^8+0.2215403975e-19*T^9-0.4589870880e-21*T^10+0.5389509965e-23*T^11, C(T) = .3000000000-0.4811931399e-1*T+0.3859529419e-2*T^2-0.2063396439e-3*T^3+0.8267055682e-5*T^4-0.2631388485e-6*T^5+0.6554216926e-8*T^6-0.5551714044e-10*T^7-0.1520371564e-10*T^8+0.2777703195e-11*T^9-0.3923000816e-12*T^10+0.4970214040e-13*T^11, DD(T) = .2100000000-.1140840000*T+0.3730565453e-1*T^2-0.8568110732e-2*T^3+0.1494601849e-2*T^4-0.2091819518e-3*T^5+0.2441368229e-4*T^6-0.2442660302e-5*T^7+0.2138541461e-6*T^8-0.1664279070e-7*T^9+0.1165691616e-8*T^10-0.7422668732e-10*T^11, E(T) = .1400000000+.1469441693*T-0.3995870043e-1*T^2+0.8709964152e-2*T^3-0.1500287681e-2*T^4+0.2093639084e-3*T^5-0.2441847491e-4*T^6+0.2442759169e-5*T^7-0.2138545263e-6*T^8+0.1664257448e-7*T^9-0.1165658906e-8*T^10+0.7422262550e-10*T^11, F(T) = .7000000000+.6324990290*T-0.1581608397e-2*T^2+0.2679492283e-5*T^3-0.6093270788e-8*T^4+0.1272029637e-9*T^5-0.4253521584e-11*T^6+0.1189725836e-12*T^7-0.2702521878e-14*T^8+0.3074973933e-16*T^9+0.3481499101e-17*T^10-0.6346661500e-18*T^11, G(T) = .4500000000-0.2249028960e-2*T+0.5983397429e-5*T^2-0.5345061622e-7*T^3+0.2810718704e-8*T^4-0.1239204116e-9*T^5+0.4250786124e-11*T^6-0.1189706297e-12*T^7+0.2702520657e-14*T^8-0.3074973865e-16*T^9-0.3481499102e-17*T^10+0.6346661500e-18*T^11}

 

 

 

 

 

 

 

 


 

Download odeRes2.mw

@Christian Wolinski 

Mine is set to [10,10] which is the default.

Until Maple2019 any matrix bigger than the interface(rtablesize) setting  wouldn't display at all (unless one reset the rtablesize). By default in Maple 2019 it now displays the matrix up to the rtablesize setting with the other entries elided, and just the little indicator at the bottom which tells you the actual matrix size - so the "15*9 Matrix" in  the output of my worksheet above

@Christian Wolinski 

One can get the output matrix directly

restart;
incog:=[theta[1, 1], theta[1, 2], theta[2, 1], theta[2, 2], theta[2, 6],
        theta[3, 0], theta[3, 3], theta[3, 4], theta[3, 5]
       ]:
eq:=[ 1, theta[1, 1]+theta[2, 2]+theta[3, 3], -theta[1, 1]-theta[2, 2],
      theta[2, 6]*theta[3, 5], -theta[1, 1]*theta[3, 3]-theta[2, 2]*theta[3, 3],
     -theta[1, 1]*theta[2, 6]*theta[3, 5]+theta[1, 2]*theta[2, 6]*theta[3, 4],
      theta[1, 1]*theta[2, 2]*theta[3, 3]-theta[1, 2]*theta[2, 1]*theta[3, 3]+theta[1, 2]*theta[2, 6]*theta[3, 0]
    ]:
myProc:= proc(P, V)
              local c, C, v;
              coeffs(P, V, 'c');
              seq([seq(`if`(type(C, dependent(v)),1,0), v = V)], C = c);
         end proc:
Matrix(map(myProc, eq, incog));

_rtable[18446744074332602366]

(1)

 


 

Download buildMat.mw

@acer 

I'm thinking about building a new workstation, using an AMD processor. (I've never gone the AMD route before - but you know - current price/performance etc!)

The fact that Maple itself appears to use the Intel MKL has now got me (slightly) worried, since (many years ago?)  this library was nicknamed "crippleAMD".

According to https://en.wikipedia.org/wiki/Math_Kernel_Library,

As of 2019, MKL, which remains the choice of many pre-compiled Mathematical applications on Windows (such as NumPy, SymPy, and MATLAB), still significantly underperforms on AMD CPUs with equivalent instruction sets

Would I really get a performance penalty if my (future) workstation uses an AMD processor??

@vv 

was actually much more interested in the discrepancy which arises between 'physics' versions

@Tegewaldt 

which command in the following you do not understand, or find "needlessly complicated"

  restart;
#
# Define a function
#
  f:= x-> 4*sin(x)+2*x^2+x+3;

proc (x) options operator, arrow; 4*sin(x)+2*x^2+x+3 end proc

(1)

#
# Return the value of this function
#
  f(p);

4*sin(p)+2*p^2+p+3

(2)

#
# Return the value of a function at some
# value of the independent variable
#
  f(2);

4*sin(2)+13

(3)

#
# Return the derivative of the function.
# Note that this derivative is itself a
# funcion
#
  D(f);

proc (x) options operator, arrow; 4*cos(x)+4*x+1 end proc

(4)

#
# Return the derivative of the funcion at a
# some value of the independent variable
#
  D(f)(2);

4*cos(2)+9

(5)

#
# Define a "new" function which contains the
# derivative of the original function
#
  g:= (x, y)-> D(f)(x)*y;

proc (x, y) options operator, arrow; (D(f))(x)*y end proc

(6)

#
# Return the value of the "new" function
#
  g(p,q);

(4*cos(p)+4*p+1)*q

(7)

#
# Return the value of the "new" function for
# a few values of the independent variables
#
  g(p,2);
  g(3,q);
  g(3,2);

8*cos(p)+8*p+2

 

(4*cos(3)+13)*q

 

8*cos(3)+26

(8)

 

Download basicDiff.mw

@Preben Alsholm 

On 64-bit Windows 7, I just updated to Physics 428.

The pop-up showed the usual progress bar, and the associated annotation reported "downloading"," installing" then "Done!".

Subsequent Physics:-Version() check in Maple gave the message

`The "Physics Updates" version in the MapleCloud is 428 and is the same as the version installed in this computer, created 2019, September 22, 7:8 hours, found in the directory C:\\Users\\TomLeslie\\maple\\toolbox\\2019\\Physics Updates\\lib\\`

as expected.

So I had no issues with this update

@minhhieuh2003 

was introducd in Maple 2015, and updated in Maple 2016. Exactly what version are you using?

@Jayaprakash J 

so you will have to specify the requirement more clearly.

Possibly??? the command

P:=unapply(doDiff(f,x,n),n);

@Melvin Brown 

This is the same error message which I obtained in a earlier post - and I suggested why this error message might occur. Try reading my earlier post

1 2 3 4 5 6 7 Last Page 2 of 147