tomleslie

Interesting problem...

Answered: tomleslie 13876

September 29 2019

3 11

One can use brute force, because it is relatively trivial to code, then check to see how bad timings get as the number of odd inputs and even addends is increased.

On my machine, the brute force approach in the attached produced the following approximate timings

#    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

Out of idle curiosity, with Nodd=10000 and Neven =10000, there were 4534 "successful" odd numbers and 463 "failures".

A brief examination of the "failures" up to these limits indicates that they always have 3 as a prime factor - coincidence???

NB It would not surprise me if there were a much more efficient way to implement this test

>

  restart;
#
# 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 f:=NumberTheory:-PrimeFactors(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)]):

memory used=1.73GiB, alloc change=40.00MiB, cpu time=14.98s, real time=15.04s, gc time=577.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) ];

(1)

>

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

(2)

>

Download prmeProb.mw

A couple of things to try...

Answered: tomleslie 13876

September 29 2019

0 1

(I'm on a Windows machine, so I can't debug/reproduce fully)

Open relevant help page with ?Physics,Geodesics
In the help page menu system, ensure that "View->Display Examples with 2D math input" is unchecked. In other words all the executable math in the help page should be displayed as 1-D input
In the help page menu system, use View->Open Page as Worksheet. which does exactly what it says. It will open a worksheet called [Physics, Geodesics], which you can execute in the standard way
Execute this worksheet one "group" at a time - What happens?
If this "works" you can try steps 1-4 gain, except at step (2) above, ensure that "View->Display Examples with 2D math input" is checked. - Again, What happens?

No restriction...

Answered: tomleslie 13876

September 27 2019

0 1

the conditions f(0)=0 and f'(0)=0 do not impose any restrictions on the values of the parameters n[1], n[2], n[3].

See the attached

>

   restart;
#
# Define the function
#
  f:=x-> 0.8680555556e-1*x*(3.262720*x^5*n[2]-.63360*x^5*n[1]^2+2.69184*x^5*n[1]-.480*x^3*n[1]^2+2.5920*x^3*n[2]+5.760*n[1]*x+1.920*n[2]*x^2+0.96e-1*x^4*n[1]^3+3.04320*x^4*n[2]-.5280*x^4*n[1]^2-0.32e-1*x^5*n[1]^4+.2976*x^5*n[1]^3-.1184*x^5*n[2]^2+.5376*x^5*n[3]-0.576e-1*x^4*n[1]+.576*x^4*n[3]-6.3360*x^3-6.96960*x^4-7.349760*x^5+11.520-.192*x^5*n[1]*n[3]-.7104*x^4*n[1]*n[2]+.2528*x^5*n[1]^2*n[2]-1.81824*x^5*n[1]*n[2]):
#
# Rewrite the above equation in a somewhat "neater"
# form by collecting powers of 'x', just to make
# it easier to read/interpret
#
  collect(expand(f(x)), x);
#
# Evaluate f(x) at x=0
#
  eval(f(x), x=0);
#
# Differentiate the above expression wrt 'x' and
# evaluate at x=0
#
  D(f)(0);

(1)

>

Download poly.mw

Like this...

Answered: tomleslie 13876

September 25 2019

3 11

which

plots all dependent variables
produces a matrix of dependent variable values for T=0..20

>

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

>

ans := dsolve([ODEs, bcs], numeric); for j in indets([ODEs], function(name)) do plots:-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 = 18446744074595705910)

(1)

>

Download odeRes.mw

Another way...

Answered: tomleslie 13876

September 25 2019

1 0

just because I like alternatives

>

  restart;
#
# Random matrix for test purposes
#
  A:=LinearAlgebra:-RandomMatrix(10):
#
# Sort columns of A in ascending order
#
  B:=Matrix([seq(sort(A[..,j]), j=1..10)], scan=rows);
#
# Define "finder" function
#
  g:=(mat, col, val)-> ListTools:-SelectFirst
                       ( x-> evalb( x>val ),
                         mat[..,col],
                         output=indices
                       ):
#
# Usage example
#
# In the matrix B, output the row index of the
# first entry in column 2 which is >50
#
  g(B, 2, 50);
#
# Returns NULL if no entry exceeds supplied value
#
  g(B, 2, 100);

(1)

>

Download searchMat.mw

One way...

Answered: tomleslie 13876

September 25 2019

1 0

is shown in the attached

Download altSum.mw

Remove the cr*p...

Answered: tomleslie 13876

September 23 2019

0 0

If I strip out all the redundant stuff you neither want nor need the I am left wiith the worksheet shown in the attached.

There are two reasons why no solution can be obtained from this worksheet, both of which you will have to fix

You have six ODES , but seven independent variables, namely A(T), B(T), C(T), G(T), R(T), S(T), I(T). Rouben has already pointed out that one of these dependent variables ( ie I(T) ) has a name which will cause problems. You need to work out if you really do have seven dependent variables: if you do, then don;'t call one of them I(T) - any other name will be fine: eg P(T), Q(T), Z(T), whatever, just not I(T). However if you do have seven dependent variables, then either
1. One of these is going to have to be defined in terms of the other six, or
2. You are going to need another equation
The variable mu__r is used in a couple of these equations and is never defined

see the attached

>

restart:

>	# # Define gamma as local (don't like doing this!) # local gamma:

>	A__h:= .18: beta__1:= 0.8354e-1: psi:= 0.7258e-1: mu__r:= 0.51e-1: sigma:= .165: alpha:= .65: xi:= 0.5e-1: gamma:= .131: A__b:= .7241: beta__o:= 0.65e-1: mu__b = 0.17e-1:

>

  odes:= diff(S(T), T) = A__h-psi*beta__1*S(T)*G(T)-mu__r*S(T),
          diff(G(T), T) = psi*beta__1*S(T)*G(T)-sigma*psi*beta__1*R(T)*B(T)-(alpha+xi+mu__r)*G(T),
          diff(A(T), T) = alpha*G(T)-(gamma+mu__r)*A(T),
          diff(R(T), T) = gamma*A(T)+sigma*psi*beta__1*R(T)*B(T)- mu__r*R(T),
          diff(C(T), T) = A__b - psi*beta__o*C(T)*I(T)- mu__b*C(T),
          diff(B(T), T) = psi*beta__o*C(T)*I(T)- mu__b*B(T):
  ics:= S(0)=100, G(0)=190, A(0)=45, R(0)=20, C(0)=35, B(0)=25:

>	# # Solve system # ans:= dsolve([ odes, ics ], numeric);

Warning, The use of global variables in numerical ODE problems is deprecated, and will be removed in a future release. Use the 'parameters' argument instead (see ?dsolve,numeric,parameters)

(1)

>

Download odeProb.mw

More like this...

Answered: tomleslie 13876

September 23 2019

1 0

maybe?

You can adjust the decay rate by adding a parameter to the exponent in the definition of the function g(). I used 0.5, just for illustration.

Download staircase2.mw

Adding a "staircase"...

Answered: tomleslie 13876

September 23 2019

1 1

It seems as if you want to add a "staircase" component to a function whcih already exists (I might be wrong!)

You might want to consider the "toy" example in the attached which adds 150 to the function 20*sin(z), at z=24, 48, etc

Download staircase.mw

Definitely something funny (again!)...

Answered: tomleslie 13876

September 23 2019

0 2

In the attached, I'm pretty sure that I am comparing the results using Physics:-Version(426) and Physics:-Version(429), although the message returned by the Physics:-Version() command is somewhat ambiguous.

In any case the two execution groups return different answers.

The first, using Physics:-Version(426) returns what I would consider to be the "correct" answer and the pdetest() command comfirms it (ie returns 0).

The second using Physics:-Version(429), returns the answer in your original post, which I agree is wrong!

Are you some kind of official beta-tester for Maplesoft??

>	restart; Physics:-Version(426); restart; Physics:-Version(); pde:=diff(u(x,t),t$2)=4diff(u(x,t),x$2); ic:=u(x,0)=0,D[2](u)(x,0)=sin(x)^2; bc:=u(-Pi,t)=0,u(Pi,t)=0; sol:=pdsolve([pde,ic,bc],u(x,t)); pdetest(sol,pde); simplify(diff(rhs(sol),t$2)-4diff(rhs(sol),x$2));

Warning, this package updates content shipped in a standard Maple install. Use the 'restart' command to clear your session before using these commands.

$Kernel(`The "Physics Updates" version "426" is installed but is not active. The active version of Physics is within the library C:\Users\TomLeslie\maple\toolbox/2019/Physics Updates/lib\Physics Updates.maple, created 2019, September 23, 10:13 hours`), [`The "Physics Updates" version "426" is installed but is not active. The active version of Physics is within the library C:\Users\TomLeslie\maple\toolbox/2019/Physics Updates/lib\Physics Updates.maple, created 2019, September 23, 10:13 hours`]$

$`The "Physics Updates" version "426" is installed but is not active. The active version of Physics is within the library C:\Users\TomLeslie\maple\toolbox/2019/Physics Updates/lib\Physics Updates.maple, created 2019, September 23, 10:13 hours`$

u(x, t) = Sum((1/2)*(Int(sin(n*x)*sin(x)^2, x = -Pi .. Pi))*sin(2*n*t)*sin(n*x)/((Int(sin(n*x)^2, x = -Pi .. Pi))*n), n = 1 .. infinity)

(1)

>	restart; Physics:-Version(429); restart; Physics:-Version(); pde:=diff(u(x,t),t$2)=4diff(u(x,t),x$2); ic:=u(x,0)=0,D[2](u)(x,0)=sin(x)^2; bc:=u(-Pi,t)=0,u(Pi,t)=0; sol:=pdsolve([pde,ic,bc],u(x,t)); pdetest(sol,pde); simplify(diff(rhs(sol),t$2)-4diff(rhs(sol),x$2));

Warning, this package updates content shipped in a standard Maple install. Use the 'restart' command to clear your session before using these commands.

$Kernel(`The "Physics Updates" version "429" is installed but is not active. The active version of Physics is within the library C:\Users\TomLeslie\maple\toolbox/2019/Physics Updates/lib\Physics Updates.maple, created 2019, September 23, 10:14 hours`), [`The "Physics Updates" version "429" is installed but is not active. The active version of Physics is within the library C:\Users\TomLeslie\maple\toolbox/2019/Physics Updates/lib\Physics Updates.maple, created 2019, September 23, 10:14 hours`]$

$`The "Physics Updates" version "429" is installed but is not active. The active version of Physics is within the library C:\Users\TomLeslie\maple\toolbox/2019/Physics Updates/lib\Physics Updates.maple, created 2019, September 23, 10:14 hours`$

(2)

>

Download WaveEq.mw

Yet another answer...

Answered: tomleslie 13876

September 22 2019

2 0

just because alternatives are good!

>

restart;
with(inttrans):
alias( U__1(s)=laplace(u__1(t), t, s),
       U__2(s)=laplace(u__2(t), t, s),
       Y__1(s)=laplace(y__1(t), t, s),
       Y__2(s)=laplace(y__2(t), t, s)
     ):
sys:=[ 2*diff(y__1(t),t)=-2*y__1(t)-3*y__2(t)+2*u__1(t),
       diff(y__2(t),t)=4*y__1(t)-6*y__2(t)+2*u__1(t)+4*u__2(t)
     ];
lsys:=eval( laplace~(sys, t,s), [y__1(0)=0, y__2(0)=0]);

(1)

>	eq1:=isolate( subs(isolate(lsys[1], Y__2(s)),lsys[2]),Y__1(s)); eq2:=simplify(isolate( subs(eq1, lsys[2]),Y__2(s)));

Y__1(s) = (-2*U__1(s)+4*U__2(s)-(2/3)*s*U__1(s))/(-(2/3)*s^2-(14/3)*s-8)

Y__2(s) = ((2*s+6)*U__1(s)+4*U__2(s)*(s+1))/(s^2+7*s+12)

(2)

>	G_11:= simplify(eval(eq1, U__2(s)=0))/U__1(s); G_12:= simplify(eval(eq1, U__1(s)=0))/U__2(s); G_21:= simplify(eval(eq2, U__2(s)=0))/U__1(s); G_22:= simplify(eval(eq2, U__1(s)=0))/U__2(s);

Y__1(s)/U__1(s) = 1/(s+4)

Y__1(s)/U__2(s) = -6/(s^2+7*s+12)

Y__2(s)/U__1(s) = 2/(s+4)

Y__2(s)/U__2(s) = 4*(s+1)/(s^2+7*s+12)

(3)

>

Download lap.mw

Errr...

Answered: tomleslie 13876

September 22 2019

1 1

Isn't this a "straightforward" z-gradient scheme with a defined zrange, as in the attached (which has the additional benefit of working in Maple 2015).

Or am I missing something??

PS Plot colours render a lot better in Maple than they do on this site!

>

  restart:
  interface(version);
  with(plots):
  with(Statistics):
  randomize():
  N := 10:
  P := 3:
  A := Sample(Uniform(0, 1), [N, P]):
  C := CorrelationMatrix(A);
  matrixplot
  ( C,
    heights=histogram,
    axes=frame,
    gap=0.25,
    colorscheme=[ "zgradient",
                  [ "Blue", "White", "Red" ],
                  zrange=-1..1
                ],
    orientation=[0, 0, 0],
    lightmodel=none,
    tickmarks=[[seq(i+1/2=A||i, i=1..P)], [seq(i+1/2=A||i, i=1..P)], default],
    labels=[("")$3]
  );

C := Matrix(3, 3, {(1, 1) = 1.0, (1, 2) = -.1563106801276543, (1, 3) = -0.6792568803622306e-1, (2, 1) = -.1563106801276543, (2, 2) = 1.0, (2, 3) = -.37067762598813236, (3, 1) = -0.6792568803622306e-1, (3, 2) = -.37067762598813236, (3, 3) = 1.0}, datatype = float[8])

>

Download zgrad.mw

Defintely something funny...

Answered: tomleslie 13876

September 21 2019

1 0

In the attached worksheet (on my machine)

it starts with Physics:-Version 419 and everything works
it sets the Physics:-Version to 426 repeats the commands, and all integrals return unevaluated!!
it sets the Physics:-Version back to 419 and everything works again

If you re-execute this worksheet, your experience will be somewhat different, because you will be starting with Physics:-Version 426, so the first two sets of commands *should* give the same "unevaluated integrals. But what happens with the last set of commands where you should be using Physicss:-Version 419??

>

restart;

>

version()

User Interface: 1399874

Kernel: 1399874
Library: 1399874

(1)

>	interface(version)

(2)

>	Physics:-Version();

$`The "Physics Updates" version "419" is installed but is not active. The active version of Physics is within the library C:\Users\TomLeslie\maple\toolbox\2019\Physics Updates\lib\Physics Updates.maple, created 2019, September 21, 9:30 hours`$

(3)

>

restart;

>	int(exp(x),x=0..1)

(4)

>	int(sin(n*x),x=0..Pi)

(5)

>	int(tan(x),x=0..Pi)

(6)

>	int(cos(x),x=0..1)

(7)

>	int(sin(x),x=0 .. Pi)

(8)

>	int(cos(x),x)

(9)

>	int(x,x=0 .. 1)

(10)

>	Physics:-Version(426);

Warning, this package updates content shipped in a standard Maple install. Use the 'restart' command to clear your session before using these commands.

$`The "Physics Updates" version "426" is installed but is not active. The active version of Physics is within the library C:\Users\TomLeslie\maple\toolbox/2019/Physics Updates/lib\Physics Updates.maple, created 2019, September 21, 9:34 hours`$

(11)

>

restart;

>	int(exp(x),x=0..1)

(12)

>	int(sin(n*x),x=0..Pi)

int(sin(n*x), x = 0 .. Pi)

(13)

>	int(tan(x),x=0..Pi)

int(tan(x), x = 0 .. Pi)

(14)

>	int(cos(x),x=0..1)

int(cos(x), x = 0 .. 1)

(15)

>	int(sin(x),x=0 .. Pi)

int(sin(x), x = 0 .. Pi)

(16)

>	int(cos(x),x)

(17)

>	int(x,x=0 .. 1)

(18)

>	Physics:-Version(419);

Warning, this package updates content shipped in a standard Maple install. Use the 'restart' command to clear your session before using these commands.

$`The "Physics Updates" version "419" is installed but is not active. The active version of Physics is within the library C:\Users\TomLeslie\maple\toolbox/2019/Physics Updates/lib\Physics Updates.maple, created 2019, September 21, 9:34 hours`$

(19)

>

restart;

>	int(exp(x),x=0..1)

(20)

>	int(sin(n*x),x=0..Pi)

(21)

>	int(tan(x),x=0..Pi)

(22)

>	int(cos(x),x=0..1)

(23)

>	int(sin(x),x=0 .. Pi)

(24)

>	int(cos(x),x)

(25)

>	int(x,x=0 .. 1)

(26)

>

Download PhysVer.mw

Or...

Answered: tomleslie 13876

September 20 2019

0 0

If you want better control over the x-values at which the function is evaluated, use the attached.

You will have to change the filename in the Export() command to something appropriate for your machine

>

  restart;
  f:= x->x^3-2*x:
  startVal:= 1;
  stopVal:= 10;
  spacing:= 0.1;
  M:= Matrix
      ( 1+round((stopVal-startVal)/spacing),
        2,
        (i,j)-> `if`( j=1,
                      startVal+(i-1)*spacing,
                      f( startVal+(i-1)*spacing)
                    )
      );
  Export("C:/Users/TomLeslie/Desktop/data.csv", M);