10 years, 7 days

## Maybe...

the attached is what you want

If not you are going to have to explain your issue more clearly

 >
 >
 >
 >
 >
 > g := x -> piecewise(0 < x and x <= 10, x, 10 < x and x <= 20, -x+20);
 (1)
 >
 >
 (2)
 >
 >
 >
 >
 >
 > # # The following plots x-M*exp(-mu*x) for x=1. # Fot other values of 'x' change the final option # to x=2,2,-1, whatever #   plots:-pointplot( eval~(x-M*~exp(-mu*x), x=1));
 >

## Like this maybe...

Achieving the attached convinced me that all calculations should be done "unitless" - and (only if strictly necessary) should units be applied at the end!

I'm not even sure that it is "complete" becuase the units of the expression returned by dsolve() ought to be 'm', but the expression *appears* to be "unitless"

Anyhow, for what it is worth

 > restart; with(Units): assume(t>0);
 > T := 60*Unit('s'); mass := 80*Unit('kg'); b__1 := 13*Unit('kg')/Unit('s'); g := -9.81*Unit('m')/Unit('s')/Unit('s'); b := piecewise( 0 <= t*Unit('s') and t*Unit('s') < T,                 b__1,                 T<= t*Unit('s') and t*Unit('s') < infinity*Unit('s'),                 7*b__1               ); Acceleration := diff(h(t)*Unit('m'), t*Unit('s'), t*Unit('s')) = g - b*diff(h(t)*Unit('m'), t*Unit('s'))/mass; IC := h(0)*Unit('m') = 4000*Unit('m'),       D(h)(0)*Unit(('m')/('s')) = 0.1*Unit('m'/'s'); solution := dsolve({IC, Acceleration}, h(t));
 (1)
 > fsolve( rhs(solution), t); plot( rhs(solution),       t = 0 .. 150,       y = 0 .. 4100,       axis[2] = [gridlines = 30],       axis[1] = [gridlines = 20],       labels = ["t [s]", "h(t) [m]"],       labelfont = ["Times", 15],       color = "red", size = [1000, 1000]     )
 >

## A "slow" animation...

as in the attached. On this site, the animation appears to be runniing at 10 frames/sec (the default rate), but if you download the worksheet, you can use the worksheet animation toolbar to

1. set the update rate to 1 frame/sec
2. decide whether you wnat the animation to play a singel cycle or continuously
3. manuallly step frame-by-frame
 > restart;   rnge:=-1..1:   fcns:=[x, x^2, 1-x]:   plots:-display( [ seq                     ( plot                       ( fcns[j], x=rnge),                       j=1..3                     )                   ],                   insequence=true                 );
 >
 >

## Interesting problem...

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

## A couple of things to try...

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

1. Open relevant help page with ?Physics,Geodesics
2. 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
3. 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
4. Execute this worksheet one "group" at a time - What happens?
5. 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...

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

## Like this...

which

1. plots all dependent variables
2. produces a matrix of dependent variable values for T=0..20
 >
 >
 >
 >
 >
 (1)
 >
 >

## Another way...

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

## One way...

is shown in the attached

 > restart; f:=k->sum((r+1)*(r+2)*(r+3)*(r+4)*F(r+4)*(k-r+1)*F(k-r+1), r = 0 .. k) =c: g:=h->isolate( h, indets(h, function)[-1]): (g @ f)(0); (g @ f)(1); (g @ f)(2);
 (1)
 >

## Remove the cr*p...

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

1. 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
2. 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);
 (1)
 >

## More like this...

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.

 > restart; g:=z->exp(-0.5*(z-24*floor(z/24))); f:=(t)->g(t)*(t-24*floor(t/24))*150; plot(f, 0..100);
 >

## Adding a "staircase"...

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

 > restart; g:=z->20*sin(z); f:=(t)->g(t)+(floor(t/24))*150; plot(f, 0..100);
 >

## Definitely something funny (again!)...

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)=4*diff(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)-4*diff(rhs(sol),x\$2));
 (1)
 > restart; Physics:-Version(429); restart; Physics:-Version(); pde:=diff(u(x,t),t\$2)=4*diff(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)-4*diff(rhs(sol),x\$2));
 (2)
 >

## Yet another answer...

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)));
 (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);
 (3)
 >

## Errr...

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