Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 165 days
Contact tomleslie

MaplePrimes Activity


These are replies submitted by tomleslie

Yes I know...

tomleslie 13876
March 21 2020
0 0

@Carl Love 

but this suggestion is backwards.

If gridlines are 'off' by default in all plotting commands, they should be 'off' by default on this site

Idle curiosity...

tomleslie 13876
March 19 2020

Why don't you just read the help which comes packaged with yout Maple installation?

All  you have to do is click the 'help' on the Maple worksheet toolbar

Then you can do interesting things like

  1. toggle the examples in the help between 1D and 2D maple input
  2. convert the help file to an executable worksheet

Precisely what do you gain by looking for something on the web, then downloading it - when you already have the same capabilities locally, just by virtue of having installed Maple?

 

 

Hmmmm...

tomleslie 13876
March 19 2020
0 0

@djc 

I still cannot reproduce.

The only difference between out Maple installations is that I am running 64-bit WIndows 7 Home Premium, while you are running 64 bit Windows 10 Enterprise. Other respondents so far

ACER gets the 'correct' results on 64-bit Linux
Preben gets the correct results on 64- bit Windows - version unspecified. It might be useful to know his Windows version!

So the anomaly at the moment would seem(??) to be that you are running Windows 10.

It would obviously be helpful if someone here running Maple 2020 on a Windows 10 system could double-check your original code - because if the issue is related to a specific Windows version, then I suspect that only MapleSoft will be able to help you!

Use view option...

tomleslie 13876
March 19 2020
0 0

@maple2015 

to zoom in on the end points of the second branch plot with the code below- you will see that the boundary conditions are met

plots:-odeplot(Sol2, gridlines= false);
plots:-odeplot(Sol2, gridlines= false, view=[19.95..20.05, -1000..0]);
plots:-odeplot(Sol2, gridlines= false, view=[20.95..21.05, -400..0]);

 

Specify OS and Maple version...

tomleslie 13876
March 18 2020
0 0

because everything seems o be working OK for me - see the attached

> 

restart;
kernelopts(version);
is(0 <= (a - b)^2 + (c - d)^2) assuming real;
is(x = 0) assuming (0 < abs(x));

`Maple 2020.0, X86 64 WINDOWS, Mar 4 2020, Build ID 1455132`

  

true

  

false

 (1)
> 

 


Download newExps.mw

 

 

The attached...

tomleslie 13876
March 18 2020
0 0

@willykz1 

Maple file will solve your problem and generate the attached xlxs file.

In order to make this function correctly you will have to change the string 'pathname' to something appropriate for your machine

> 

restart
> 

g := 9.81; m := 0.27e-2; h := 0.73e-3; a := 0.18e-4; xi := 0; yi := 0; zi := 0.5e-1; vxi := 2.444; vyi := 0; vzi := 3.244; wx := 0; wy := 0; wz := 0
> 

eqx := m*(diff(x(t), t, t)) = -h*sqrt((diff(x(t), t))^2+(diff(y(t), t))^2+(diff(z(t), t))^2)*(diff(x(t), t))+a*(wy*(diff(z(t), t))-wz*(diff(y(t), t))); eqy := m*(diff(y(t), t, t)) = -h*sqrt((diff(x(t), t))^2+(diff(y(t), t))^2+(diff(z(t), t))^2)*(diff(y(t), t))+a*(wz*(diff(x(t), t))-wx*(diff(z(t), t))); eqz := m*(diff(z(t), t, t)) = -m*g-h*sqrt((diff(x(t), t))^2+(diff(y(t), t))^2+(diff(z(t), t))^2)*(diff(z(t), t))+a*(wx*(diff(y(t), t))-wy*(diff(x(t), t))); ei := x(0) = xi, y(0) = yi, z(0) = zi, (D(x))(0) = vxi, (D(y))(0) = vyi, (D(z))(0) = vzi
> 

F := dsolve({ei, eqx, eqy, eqz}, {x(t), y(t), z(t)}, numeric)
> 

tEnd := 10; plots:-odeplot(F, [x(t), y(t), z(t)], t = 0 .. tEnd)

 
> 

tStep := .1; pathname := "C:/Users/TomLeslie/Desktop/"; fname := "aTest.xlsx"; M := Matrix([`~`[convert](`~`[lhs](F(0)), string), seq(`~`[rhs](F(j)), j = 0 .. tEnd, tStep)]); ExcelTools:-Export(M, cat(pathname, fname), "all"); ExcelTools:-Export(`<|>`(M[() .. (), 1], M[() .. (), 2], M[() .. (), 4], M[() .. (), 6]), cat(pathname, fname), "partial")
> 

``

Download odeToXL.mw

aTest.xlsx

It isn't too difficult...

tomleslie 13876
March 16 2020

@Axel Vogt 

to scale everything so that N=1 - see the attached.

However the fundamental problem with this scaling is how one then defines the "classic" inital conditions - ie that there is 1 infected person and everyone else in the population is susceptible.

Having defined N=1 in the attached, I further defined the initial infection as I(0)=1e-08, and the initial susceptibles as S(0)=1-1e-08: but (by implication) this is actualy defining the total population as 1e-08

NB I haven;t updated the comments in the attached, so wherever it says something like "the number of people who are susceptible"  you have to interpret as "fraction of the population who are susceptible"

> 

  restart;
  with(plots):
#
# Change the imaginaary unit, just because
# most models use I() to designate the number
# of the population who are infected and I
# decided to stick with this convention
#
  interface(imaginaryunit=J):
#
# Population
#
  N:=1:
> 

####################################################
# S(t) represents the number of people who are
# susceptible to the disease
#
# I(t) represents the number of people who are
# infected
#
# R(t) represents the number of people who have
# reached a resolution - in other words they have
# either recovered or are dead!. With either of
# these resolutions, they are no longer infected or
# susceptible! If 'x' is the death rate, then
#
#      x*R(t) people are dead
#  (1-x)*R(t) people have recovered
#
# For Covid-19 the death rate is etimated anywhere
# between 1% and 5% - so 3% is probably about as good
# a guess as you are going to get
#
####################################################
#
# N represents the total population whihc is assumed
# constant - and hence neglects "normal" birth and
# death (unrelated to the epidemic) rates
#
# myBeta is the infection rate per unit time - in
# other words how many people per day will get
# the disease from a single infected person. Note that
# this can be reduced by "social isolation". If the
# infected person doesn't get to interact with the
# susceptible population, then the number of people
# who can be infected goes down
#
# myGamma is the "recovery" rate per unit time.
# Another way to look at this is that for 1/myGamma
# units of time an individual has the capability
# to infect others. After that time the infected
# individual has either recoverd, or is dead. Either
# way they no longer infect anyone else.
#
  odes:= diff(S(t),t) = -myBeta*I(t)*S(t)/N,
         diff(I(t),t) = myBeta*I(t)*S(t)/N-myGamma*I(t),
         R(t)=N-S(t)-I(t):

#
# At time=0, there has to be (at least) one person
# infected (otherwise we can't start an epidemic),
# and hence N-1 people are susceptible
#
  ics:= S(0)=N-1.0e-08, I(0)=1.0e-08:
#
# Solve the system
#
  sol:=dsolve( [odes, ics],
               parameters=[myBeta, myGamma],
               numeric
             ):
> 

###################################################################
# The first set of examples illustrate that for a given ratio of
# infection_rate/recovery_rate (I have used 2 for this ratio) the
# only thing which changes in the evolution of the epidemic is the
# timescale - other than timescale the numbers stay the same
#
# Notice that "herd immunity" kicks in. For this comibination of
# numbers, about 80% of the population are going to get it. 20% of
# the susceptible population never do.
#
# NB note the different timescales in these plots
#
###################################################################
#
# A high infection rate with a short period during which
# infection can take place
#
  sol(parameters=[1.0, 0.5]):
  odeplot( sol,
           [ [t, S(t)], [ t, I(t)], [t, R(t)] ],
           t=0..100
         );
#
# A moderate infection rate with a moderate period during which
# infection can take place
#
  sol(parameters=[ 0.5, 0.25]):
  odeplot( sol,
           [ [t, S(t)], [ t, I(t)], [t, R(t)] ],
           t=0..200
         );
#
# A low infection rate with a long period during which
# infection can take place
#
  sol(parameters=[0.1, 0.05]):
  odeplot( sol,
           [ [t, S(t)], [ t, I(t)], [t, R(t)] ],
           t=0..1000
         );

 

 

 
> 

###################################################################
#
# A more realistic(?) scenario is where the recovery rate is fixed.
# You get the disease: you are asymptomatic (but infectious at some
# rate) for a period. Then symptomatic (and still infectious) for
# another period. This is a function of the infection, and something
# over whihc no-one has any control.
#
# However, the infection rate varies substantially between these two
# periods. In the second phase you have probably retired to bed and
# are infecting no-one. You probably infected your "carers" (ie
# family) during the first period and can't infect them again!
#
# Estimates for covid-19 suggest that you will be infectious, but
# asymptomatic, for about 5 days; then infectious and symptomatic
# for about another 5 days. But for this second phase you (probably)
# won't induce any new infections. (You've already infected close
# family in phase one). This means that the best guess on the total
# time you are infectious is (probably?) ~5 days
#
# So an obvious thing to look at is a recovery rate of 0.2 (ie the
# reciprocal of the time for which one is infectious) and then play
# with the infection rate numbers - and boy do these make a big
# difference!
#
# A high infection rate
#
  sol(parameters=[1.0, 0.2]):
  odeplot( sol,
           [ [t, S(t)], [ t, I(t)], [t, R(t)] ],
           t=0..100
         );
#
# A moderate infection rate
#
  sol(parameters=[0.5, 0.2]):
  odeplot( sol,
           [ [t, S(t)], [ t, I(t)], [t, R(t)] ],
           t=0..200
         );
#
# A low infection rate
#
  sol(parameters=[0.25, 0.2]):
  odeplot( sol,  
           [ [t, S(t)], [ t, I(t)], [t, R(t)] ],
           t=0..1000
         );

 

 

 
> 

 

Download covid2.mw

My $0.02...

tomleslie 13876
March 15 2020

A couple of weeks ago, I put together a really simple model (based on some Wikipedia page, can't remember where)

This, with a few "test cases" is shown in the attached

This model is "simple", and many "more sophistcated" models are available. However, all of these models (including the one in the attached) are very sensitive to supplied parameters. Searching for ever more accurate sets of parameters, seems pointless.

The primary objective is not to model the epidemic more accurately.

The primary objective to is reduce the impact

The only thing that playing around with parameters in the attached actually told me is that "the powers that be" have only one adjustable knob - the infection rate - just get it as low as possible. So good hygiene, social isolation, etc, etc.

For a given resoluton (ie recovery or death) rate, which is entirely controlled by the virus, overall outcomes vary substantially based purely on the infection rate - and that is the only button we can push!

Anyhow FWIW if you feel like playing with numbers, check out the attached

> 

  restart;
  with(plots):
#
# Change the imaginaary unit, just because
# most models use I() to designate the number
# of the population who are infected and I
# decided to stick with this convention
#
  interface(imaginaryunit=J):
> 

####################################################
# S(t) represents the number of people who are
# susceptible to the disease
#
# I(t) represents the number of people who are
# infected
#
# R(t) represents the number of people who have
# reached a resolution - in other words they have
# either recovered or are dead!. With either of
# these resolutions, they are no longer infected or
# susceptible! If 'x' is the death rate, then
#
#      x*R(t) people are dead
#  (1-x)*R(t) people have recovered
#
# For Covid-19 the death rate is etimated anywhere
# between 1% and 5% - so 3% is probably about as good
# a guess as you are going to get
#
####################################################
#
# N represents the total population whihc is assumed
# constant - and hence neglects "normal" birth and
# death (unrelated to the epidemic) rates
#
# myBeta is the infection rate per unit time - in
# other words how many people per day will get
# the disease from a single infected person. Note that
# this can be reduced by "social isolation". If the
# infected person doesn't get to interact with the
# susceptible population, then the number of people
# who can be infected goes down
#
# myGamma is the "recovery" rate per unit time.
# Another way to look at this is that for 1/myGamma
# units of time an individual has the capability
# to infect others. After that time the infected
# individual has either recoverd, or is dead. Either
# way they no longer infect anyone else.
#
  odes:= diff(S(t),t) = -myBeta*I(t)*S(t)/N,
         diff(I(t),t) = myBeta*I(t)*S(t)/N-myGamma*I(t),
         R(t)=N-S(t)-I(t):

#
# At time=0, there has to be (at least) one person
# infected (otherwise we can't start an epidemic),
# and hence N-1 people are susceptible
#
  ics:= S(0)=N-1, I(0)=1:
#
# Solve the system
#
  sol:=dsolve( [odes, ics],
               parameters=[N, myBeta, myGamma],
               numeric
             ):
> 

###################################################################
# The first set of examples illustrate that for a given ratio of
# infection_rate/recovery_rate (I have used 2 for this ratio) the
# only thing which changes in the evolution of the epidemic is the
# timescale - other than timescale the numbers stay the same
#
# Notice that "herd immunity" kicks in. For this comibination of
# numbers, about 80% of the population are going to get it. 20% of
# the susceptible population never do.
#
# NB note the different timescales in these plots
#
###################################################################
#
# A high infection rate with a short period during which
# infection can take place
#
  sol(parameters=[100000000, 1.0, 0.5]):
  odeplot( sol,
           [ [t, S(t)], [ t, I(t)], [t, R(t)] ],
           t=0..100
         );
#
# A moderate infection rate with a moderate period during which
# infection can take place
#
  sol(parameters=[100000000, 0.5, 0.25]):
  odeplot( sol,
           [ [t, S(t)], [ t, I(t)], [t, R(t)] ],
           t=0..200
         );
#
# A low infection rate with a long period during which
# infection can take place
#
  sol(parameters=[100000000, 0.1, 0.05]):
  odeplot( sol,
           [ [t, S(t)], [ t, I(t)], [t, R(t)] ],
           t=0..1000
         );

 

 

 
> 

###################################################################
#
# A more realistic(?) scenario is where the recovery rate is fixed.
# You get the disease: you are asymptomatic (but infectious at some
# rate) for a period. Then symptomatic (and still infectious) for
# another period. This is a function of the infection, and something
# over whihc no-one has any control.
#
# However, the infection rate varies substantially between these two
# periods. In the second phase you have probably retired to bed and
# are infecting no-one. You probably infected your "carers" (ie
# family) during the first period and can't infect them again!
#
# Estimates for covid-19 suggest that you will be infectious, but
# asymptomatic, for about 5 days; then infectious and symptomatic
# for about another 5 days. But for this second phase you (probably)
# won't induce any new infections. (You've already infected close
# family in phase one). This means that the best guess on the total
# time you are infectious is (probably?) ~5 days
#
# So an obvious thing to look at is a recovery rate of 0.2 (ie the
# reciprocal of the time for which one is infectious) and then play
# with the infection rate numbers - and boy do these make a big
# difference!
#
# A high infection rate
#
  sol(parameters=[100000000, 1.0, 0.2]):
  odeplot( sol,
           [ [t, S(t)], [ t, I(t)], [t, R(t)] ],
           t=0..100
         );
#
# A moderate infection rate
#
  sol(parameters=[100000000, 0.5, 0.2]):
  odeplot( sol,
           [ [t, S(t)], [ t, I(t)], [t, R(t)] ],
           t=0..200
         );
#
# A low infection rate
#
  sol(parameters=[100000000, 0.25, 0.2]):
  odeplot( sol,  
           [ [t, S(t)], [ t, I(t)], [t, R(t)] ],
           t=0..1000
         );

 

 

 
> 

 

Download covid.mw

 

 

Don't worry about it...

tomleslie 13876
March 13 2020
0 0

@Christopher2222 

Acer's solution is much better: one can even plottools:-rotate() the animation if desired, which (slightly) surprised me  not sure why

See the attached

> 

  restart;
  with(plots):
  with(plottools):
  xx := [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]:
  yy := [3, 3, 3, 4, 5, 6, 8, 7, 6, 5]:
  s := CurveFitting:-Spline(xx, yy, v):
> 

  a1:= animate
       ( plot3d,
         [ [s, th, v],
           v=0..y,
           th=0..2*Pi,
           coords=cylindrical,
           style=surface
         ],
         y=0..9
       ):
  display(a1, insequence=true);

 
> 

  rotate
  ( a1,
    0, Pi/2, 0
  );

 
> 

 


 

Download animsurf4.mw

 

See the help...

tomleslie 13876
March 13 2020
0 0

on ?CodeTools:-Usage

I recommend...

tomleslie 13876
March 13 2020
0 0

@Christopher2222 

tht you break the problem down into simple stages until you get "comfortable" with thought process. For example

  1. Produce a simple 2D curve
  2. Produce a simple 2-D animation of the curve in (1) above
  3. Produce a 3D animation by generating a surface of revolution, based on rotating the curve in (2) above, about the x-axis. The axis of cylindrical symmetry will be the x-axis, which makes this 3D surface simple to compare with the "profile" generated in (1) above
  4. Produce a further animation (if necessary/desired) by rotating the surface generated in (3) above by Pi/2 about the y-axis. This will place the axis of cylindrical symmetry along the z-axis

Check the attached where I have illustrated the above process in four stages.

NB this site is still rotating (some of) the animations which makes stage 3 of the above process look a bit wrong. I recommend downloading/executing the attached worksheet

> 

  restart;
  with(CurveFitting):
  with(plots):
  with(plottools):
  with(Student[Calculus1]):
> 

#
# OP's required curve
#
  x := [0, 1, 2, 3, 4, 5, 6, 7, 8, 9];
  y := [3, 3, 3, 4, 5, 6, 8, 7, 6, 5];
  s := Spline(x, y, v);plot(s, v=0..max(x)):
  plot( s, v=0..max(x) );

x := [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

  

y := [3, 3, 3, 4, 5, 6, 8, 7, 6, 5]

  

piecewise(v < 1, -(562/8109)*v^3+(562/8109)*v+3, v < 2, (2810/8109)*v^3-(1124/901)*v^2+(10678/8109)*v+6985/2703, v < 3, -(2569/8109)*v^3+(2462/901)*v^2-(53870/8109)*v+21329/2703, v < 4, -(643/8109)*v^3+(536/901)*v^2-(1868/8109)*v+235/159, v < 5, (97/153)*v^3-(7176/901)*v^2+(275764/8109)*v-39799/901, v < 6, -(11812/8109)*v^3+(21079/901)*v^2-(18787/153)*v+586978/2703, v < 7, (9671/8109)*v^3-(21887/901)*v^2+(77909/477)*v-959798/2703, v < 8, -(2545/8109)*v^3+(6617/901)*v^2-(471299/8109)*v+436898/2703, (509/8109)*v^3-(1527/901)*v^2+(115069/8109)*v-28106/901)

  

 
> 

#
# Produce a simple 2-D animation of the curve
#
   display
   ( seq
     ( plot
       ( s,
         v=0..i*max(x)/60
       ),
       i=0..60
     ),
     insequence=true
   );

 
> 

#
# Produce a "surface of revolution" by rotating the above
# about the x-axis: this means that the axis of cylindrical
# symmetry is the x-axis
#
   display
   ( seq
     ( SurfaceOfRevolution
       ( s,
         v=0..i*max(x)/60,
         output=plot,
         caption=""
       ),
       i=0..60
     ),
     insequence=true
   );

 
> 

#
# Rotate the above plot so that the axis of cylindrical
# symmetry becomes the z-axis
#
   display
   ( seq
     ( rotate
       ( SurfaceOfRevolution
         ( s,
           v=0..i*max(x)/60,
           output=plot,
           caption=""
         ),
         0, Pi/2, 0
       ),
       i=0..60
     ),
     insequence=true
   );

 
> 

 


 

Download animSurf.mws
 

I would advise against overcomplicating...

tomleslie 13876
March 12 2020
0 0

@Christopher2222 

If the plots that you desire have "cylindrical symmetry", then I would recommend that you

  1. generate a curve which has the "profile" you want
  2. use the SurfaceOfRevoution() command to generate the 3-D shape that you want
  3. use the plottools:-rotate/reflect/translate commands to "orient" the 3-D shape produced at (2) above

For example the output of the following code is quite "cute", and could probabyl(?) be generated in other ways - but these would take a lot of thinking about!

restart:
  with(plots):
  with(plottools):
  with(Student[Calculus1]):display( seq
           ( rotate
             ( SurfaceOfRevolution
               ( exp(-0.1*x)*(2+sin(x)),
                 x=0..j*Pi/10,
                 output=plot
               ),
               0, Pi/2, 0
             ),
             j=1..60
           ),
           insequence=true
          );

 

I'd be using indexed variables!...

tomleslie 13876
March 10 2020
0 0

@AndreasBuus 

and then write something like

add( n[i]*(X[i]-x)^2, i=1..4);

 

Oh dear - 2D versus 1D again...

tomleslie 13876
March 06 2020
0 0

@umbli 

This could start an argument in which I have little or no interest - but make the following observations

  1. I never use 2D input for my own work - I could give you a list of reasons, but I won't.
  2. The only time I play with 2D input is when I am fixing problems on this site for those who have problems with it.
  3. If you insist on using 2D math input then I  would suggest that if you want to multiply anything in Maple then typing the '*' character on you keyboard between the multiplicands is a good idea.

 

A compromise...

tomleslie 13876
March 03 2020
0 0

@jum 

With the example you provide, ie

sys := [.8164965809*c[1]+101.3127271*c[3]-2.479560632 = 0., .7071067810*c[1]+1.414213562*c[2]+146.5538537*c[3]-2.772453851 = 0., 4/3*c[1]+c[2] = (1/24)*sqrt(Pi)*(19+16*sqrt(2))]; ics := c[1]-2*c[2]+3*c[3] = sqrt(Pi), 4*c[2]-16*c[3] = (1/2)*sqrt(Pi), c[1]+2*c[2]+3*c[3] = sqrt(2)*sqrt(Pi); 4/3*c[1]+c[2] = (1/24)*sqrt(Pi)*(19+16*sqrt(2)); fsolve(sys);
{c[1] = 2.187662474, c[2] = 0.1573948503, c[3] = 0.00684357949}

  1. You can solve exactly for 'sys' but the solution will not be valid for 'ics'
  2. You can solve exactly for 'ics' but the solution will not be valid for 'sys'
  3. You can solve "approximately" for both 'sys' and 'ics', minmising the residual errors across all 6 equations

All three of the above options are addressed in the attached

> 

  restart:

  sys := [.8164965809*c[1]+101.3127271*c[3]-2.479560632 = 0.,
          .7071067810*c[1]+1.414213562*c[2]+146.5538537*c[3]-2.772453851 = 0.,
          4/3*c[1]+c[2] = (1/24)*sqrt(Pi)*(19+16*sqrt(2))
         ]:

  ics := [ c[1]-2*c[2]+3*c[3] = sqrt(Pi),
           4*c[2]-16*c[3] = (1/2)*sqrt(Pi),
           c[1]+2*c[2]+3*c[3] = sqrt(2)*sqrt(Pi)
         ]:
#
# Solve 'sys' exactly for c[1], c[2], c[3]
#
  ans1:=fsolve(sys);
#
# Solve 'ics' exactly for c[1], c[2], c[3]
#
  ans2:=fsolve(ics);
#
# Note that the above two solution are different so that
# one cannot simultaneously satisfy both 'sys' and 'ics'
#
# Find the best(!?!) values for c[1], c[2], c[3] which
# simultaneously minimise residual errors in both 'ics'
# and 'sys'
#
  ans3:=Optimization:-Minimize( add( (rhs(j)-lhs(j))^2,
                                      j in [sys[],ics[]]
                                   )
                              );
#
# Obtain the residual errors in each of the equations
# with the above "compromise" solution
#
  evalf(eval(rhs~([sys[],ics[]])-lhs~([sys[],ics[]]), ans3[2]));

 

{c[1] = 2.187662474, c[2] = .1573948503, c[3] = 0.684357949e-2}

  

{c[1] = 2.168050908, c[2] = .1835436051, c[3] = -0.950328156e-2}

  

[0.276259794881e-1, [c[1] = HFloat(2.1351303326338247), c[2] = HFloat(0.22489317356760152), c[3] = HFloat(0.006661317985266363)]]

  

[HFloat(0.06135772446098953), HFloat(-0.031600083046590566), HFloat(0.0025445309206322264), HFloat(0.06712591154557934), HFloat(0.09323531899385573), HFloat(-0.09827235972482651)]

 (1)
> 

 

Download comprom.mw

 

 

 

First 54 55 56 57 58 59 60 Last Page 56 of 207