tomleslie

13876 Reputation

20 Badges

15 years, 181 days

MaplePrimes Activity


These are replies submitted by tomleslie

@erik10 

Dice1() and Dice2() in a way which ensures that each time you call (say) Dice1() you are going to get a 'new' value, what would you honestly expect from an expression which contained more than one instance of (say) Dice1()?

I mean would you honestly expect the effect of

Dice1()+DIce1()

to be the same as

2*Dice1()

Get real.

Once you figure out how to address this issue, your latest problem can be solved by the code in the attached

 restart;
 with(plots):
 with(Statistics):
 randomize():

 nFaces:= 6:
 nTrials:= 10000:
 Dice1:= rand(1 .. nFaces):
 Dice2:= rand(1 .. nFaces):
 f:= ()-> Dice1()+Dice2():
 s:= [seq(f(), i = 1 .. nTrials)]:
 statplot:= z-> Histogram
                ( z,
                  binbounds = [ seq
                                ( j+1/2,
                                  j = min(z)-1 .. max(z)
                                )
                              ]
                ):
 p1:= pointplot
      ( [ seq
          ( [j, numboccur(s, j)/nTrials],
            j = min(s)..max(s)
          )
        ],
        symbol=solidcircle,
        symbolsize=14,
        color=red,
        size=[1000, 500]
      ):
 display([p1, statplot(s)]);

 

g1 := (a, b)-> piecewise( a = 1 or b = 1,
                          -4,
                          abs(b-a)
                        ):
sq := [ seq
        ( g1(Dice1(), Dice2()),
          i = 1 .. nTrials
        )
      ]:
p1:= pointplot
      ( [ seq
          ( [j, numboccur(sq, j)/nTrials],
            j = min(sq)..max(sq)
          )
        ],
        symbol=solidcircle,
        symbolsize=14,
        color=red,
        size=[1000, 500]
      ):
display( [p1,statplot(sq)]);

 

 

``

NULL


 

Download diceProb4.mw


 

 

@Josolumoh 

NOT pictures of code. THe only way I can do anything with your previous post is to retype the oced in it. Do you really think I can be bothered?

Mapleprimes gives you a way of uploading code in a manner which means everyone here can (more or less instantaneously) read, execute, edit your code.

It is the big green up-arrow in the Mapleprimes toolbar.

If you want answers, I suggest you learn to use it.

@sand15 

that I really understand the problem. Does the attached achieve your objective.

(I have changed the range of values on the two Dice to be non-overlapping so that you have a clearer idea of what is going on)

restart:
randomize():         
with(Statistics):         
Dice1 := RandomVariable(DiscreteUniform(1, 6));
Dice2 := RandomVariable(DiscreteUniform(7, 12));
zip((x,y)->[x,y], Sample~([Dice1, Dice2], 10)[]);

Dice1 := _R

 

Dice2 := _R0

 

Vector[row](%id = 18446744074405692110)

(1)

 


 

Download twoDice.mw

 

@Carl Love 

The OP's statement

How many times ( N) do we have to roll the die to have probability at least one half that all six numbers appear in our sequence

requires the median, not the mean.

Everything I have suggested previously refers to the mean number of throws, and is therefore incorrect.

The simulation part of my original worksheet(s) can be trivially modified to yield the median, by the simple expedient of "sorting"  the results of the trials, and "picking" the one in the middle. Strictly speaking the "one in the middle" only exists for an odd number of trials, so I have wimped out and given the "one in the middle" when the number of trials is odd, and the two in the middle when the number of trials is even (see the revised attached)

For a conventional 6-face die, the median value is 13 (as given by Kitonum).

I will also withdraw my comment about this being an instance of the "Coupon Collector Problem" - https://en.wikipedia.org/wiki/Coupon_collector%27s_problem, since this appears to compute the mean (not the median), rounded up to the next integer

  restart;
#
# Set number of faces on die
#
  nFaces:=6:
#
# Using logical thought
#
  sum(nFaces/k, k=1..nFaces);

147/10

(1)

#####################################################
# Using a simulation to roll the dice and counting
# how many throws have to be made before all 6 face
# values have appeared.
#
# Repeat this simulation for nTrials
#
# For a die with 6 faces nTrials~10000 seems to be
# a reasonable tradeoff between "accuracy" and
# execution time)
#
  nTrials:= 10000:
  randomize():
  roll:= rand(1..nFaces):
  nthrows:= NULL:
  for k from 1 by 1 to nTrials do
      throws:= {};
      for j from 1 do
        #
        # Keep throwing the dice. Stop when
        # all faces have occurred
        #
          throws:= throws union {roll()};
          if   numelems(throws)=nFaces
          then break
          fi
      od:
    #
    # Record how many throws it took on this trial
    # to get all nFaces
    #
      nthrows:= nthrows, j;
  od:
#
# Average number of throws before all n faces have
# occurred
#
  evalf(add([nthrows])/numelems([nthrows]));
#
# Maximum/minimum number of throws before all n
# faces have occurred. The minimum is obviously
# nFaces (if one gets *really* lucky). The
# maximum is theoretically infinite, (although
# the chance of this happening is vanishingly small!)
#
  max([nthrows]);
  min([nthrows]);
#
# Instead of computing the *mean* number of throws
# from this simulation, compute the *median* number
# of throws by sorting the data and pickig the one
# in the middle (and please don't make me do this
# for an odd number of throws!)
#
  L:=sort([nthrows]):
  if   type( numelems(L), odd)
  then L[ceil(numelems(L)/2)];
  else L[numelems(L)/2..numelems(L)/2+1]
  fi;

14.67920000

 

71

 

6

 

[13, 13]

(2)

 

 

Download diceProb3.mw

@Joe Riel 

but since I'm bad at this kind of probability, I may well be wrong!

I'm pretty sure the the "simulation" part of my worksheet is correct. Each "trial" in this simulation is to keep rolling the die until all the faces have occurred once. Having done this for a number of trials, compute the average number of rolls it takes to get all of the faces at least once. No reason why this should be an integer!

For a die with 6 faces, both the formula and the simulation given 14.7. If you feel that this number *should* be an integer, then my interpretation would be that with 14 rolls the probability of getting all 6 faces is slightly less than 1/2, and with 15 rolls, the probability of getting all 6 faces is slightly more than 1/2

Having done a little research (which I probably should have performed, before my first post), this question appears(?) to be an instance of the "Coupon Collector Problem" - see https://en.wikipedia.org/wiki/Coupon_collector%27s_problem

The attached is a slight revision of my original, so that it is easier to run for any number of faces and any number of trials in the simulation. Note that the minimum number of rolls is always equal to the number of faces on the die: the maximum number of rolls is (potentially) infinite - so theoretically my code could run "forever", although (luckily?!) the odds of this happening is very small, so I'm not trapping it

  restart;
#
# Set number of faces on die
#
  nFaces:=6:
#
# Using logical thought
#
  sum(nFaces/k, k=1..nFaces);

147/10

(1)

#####################################################
# Using a simulation to roll the dice and counting
# how many throws have to be made before all 6 face
# values have appeared.
#
# Repeat this simulation for nTrials
#
# For a die with 6 faces nTrials~10000 seems to be
# a reasonable tradeoff between "accuracy" and
# execution time)
#
  nTrials:= 10000:
  randomize():
  roll:= rand(1..nFaces):
  nthrows:= NULL:
  for k from 1 by 1 to nTrials do
      throws:= {};
      for j from 1 do
        #
        # Keep throwing the dice. Stop when
        # all faces have occurred
        #
          throws:= throws union {roll()};
          if   numelems(throws)=nFaces
          then break
          fi
      od:
    #
    # Record how many throws it took on this trial
    # to get all nFaces
    #
      nthrows:= nthrows, j;
  od:
#
# Average number of throws before all n faces have
# occurred
#
  evalf(add([nthrows])/numelems([nthrows]));
#
# Maximum/minimum number of throws before all n
# faces have occurred. The minimum is obviously
# nFaces (if one gets *really* lucky). The
# maximum is theoretically infinite, (although
# the chance of this happening is vanishingly small!)
#
  max([nthrows]);
  min([nthrows]);

14.66490000

 

63

 

6

(2)

 

 


 

Download diceProb2.mw

 

@Carl Love 
At no time, nowhere have I claimed to compute the L^2-norm

OP's first requst was for the maximum absolute error between an "exact" and an "approximate" representation of a function, for whihc I provided a solution

Op's second request was

  1. for the maximum absolute error between an "exact" and an "approximate" representation of a different function for which I provided (more-or-less) the same solution
  2.   and , in addition, "Also , can you tell me how we can find L^(2) error for same problem" -, which I completely overlooked,- my bad.

@Nia Dutta 

Fro problems such as this, there are two possible routes to a solution

  1. Find a strictly numerical solution using strictly numerical methods - for which the best approach is probably(?) to use routines from the Optimization package. Your options here are probably NLPSolve() with the "maximise=true" option, or possibly the Maximize() function. The distinction beween these two is that the latter will only work when the expression to be "maximised" is twice differentiable: the fomer will (sometimes) work when the expression to be maximised is not differentiable at all. As a general rule using an expression containing abs() is pretty much guaranteed not to produce a differentiable function if the original expression passes through zero somewhere in its domain. Thus any method which relies on "differentiability" is may(?) fail when you deliberatedly supply a function (eg abs()) which is guaranteed to be non-differentiable under certain conditions
  2. Find a "symbolic" solution using the maximize() function. Note the lower-case maximize() which is a different animal from the Optimization:-Maximize() command; the latter invokes a numerical calculation whilst the former invokes a symbolic computation. However when using the (symbolic) maximize() the issue of differentiability arises again. HAndling differentiability in the presence of an abs() function gets difficult. Under many (but not all) circumstances one can sonsider replaicng abs(x) with x^2. Now this doesn't always work, and I can think of cases where might fail. However it generally does not introduce discontinuities. Consider the "toy" example of the function x-1. This function is twice differentiable for any value of 'x': however abs(x-1) is not differentiable at x=1; but (x-1)^2 is differentiable everywhere. Thus any method whihc relies on differentiability for a solution is likely to have more trouble with an expression containing abs()
  3. I'm a pragmatist - I tend to use what works

I did it last time (more or less)

  restart;

#
# Reset digits because is the OP is using coefficients
# of about 16 digits
#
  Digits:=20;
  exact := x^4-x^3;
  app := -2.220446049250313*10^(-16)+.19676728660169784*x+4.6491983150099205*x^2
         -3.9010072748158082*x^3+2.3619056303875987*x^4-.33198038154796344*x^5;
#
# Plot the function, just for illustration
#
  plot(abs(app-exact), x = 0 .. 2);
#
# Numeric answer from optimization package
#
  numans:=Optimization:-NLPSolve(abs(app-exact), x = 0 .. 2, maximize = true);

20

 

x^4-x^3

 

-0.22204460492503130000e-15+.19676728660169784*x+4.6491983150099205*x^2-3.9010072748158082*x^3+2.3619056303875987*x^4-.33198038154796344*x^5

 

 

[6.949387511383360978, [x = 2.]]

(1)

#
# 'Symbolic' answer from maximize command (don't use abs()
# square the desired function and take the square root of
# the solution
#
  symbans:=[ maximize( (app-exact)^2, location=true, x=0..2) ][2][];
  symbans:=[symbans[1], sqrt(symbans[2])];

[{x = 2.}, 48.293986783371023147]

 

[{x = 2.}, 6.9493875113833609809]

(2)

 

NULL


 

Download getError.mw

@Josolumoh 

distributions for a couple of random variables, but you need to generate samples from these distributions in order to export them. As in

   restart;
   N:= 50:
#
# Define and sample a couple of distributions
#
   X_1:= Statistics:-Sample(Statistics:-RandomVariable(Bernoulli(1/2)),N);
   X_2:= Statistics:-Sample(Statistics:-RandomVariable(Normal(0,1)), N);
#
# Output both sample vectors as csv
#
   ExportVector( "C:/Users/TomLeslie/Desktop/testvecX1.csv", X_1, target=csv);
   ExportVector( "C:/Users/TomLeslie/Desktop/testvecX2.csv", X_2, target=csv);

Change filepaths in the above to something appropriate for your machine

@Josolumoh 

Apparently, contradicting all your earlier posts,  you no longer wish to "Export to Excel", and you wish to  "Export to csv".

Exporting a vector to csv format is  trivial.

Use the attached with your own definition of the vector V' and a filepath appropriate for your machine

restart;
V:=LinearAlgebra:-RandomVector(10);
ExportVector( "C:/Users/TomLeslie/Desktop/testvec.csv", V, target=csv);

@Josolumoh 

"Exporting vector to excel as csv"

Now I understand

Export to Excel, and

Export as csv

but why would you want to export as csv, and then (presumably) import it to Excel (ie convert from csv to xls(x), This makes no sense!

 

@acer 

I'm impressed by these last two posts for a number of reasons

  1. Identifying that (nearly?) all of the numerical coefficients in these equations are floating point approximations to rationals. Obviously if I examine carefully I can identify the occurrence of sqrt(2), but being able to identify 9.797958972=4*sqrt(6) and 37.94733192=12*sqrt(10) etc. What made you suspect?
  2. I can see the "system" which you adopted for supplying range limits on variables: ie if the second index on a variable name is '0', then impose the range -12..12; for all other variables, no range constraint is applied. But on what basis did you come up with this "system"?
  3. BTW, whilst playing with your approach, I (more-or-less) inadvertently imposed the range constraint-12..12 on all variables. fsolve() produced an answer, albeit completely different from the one you reported. That is

fsolve( [ MyEqs ],
            { seq(j = 0, j in indets([MyEqs])) },
            { seq(j = -12 .. 12, j in indets([MyEqs])) }
         );
eval([MyEqs], sol);

returns

{a[1][0] = .5050762724, a[1][1] = -8.898452132*10^(-13), a[1][2] = 4.994293575*10^(-14), a[2][0] = 0., a[2][1] = 5.029888588*10^(-13), a[2][2] = -4.994293365*10^(-14), b[1][0] = 0., b[1][1] = 0., b[1][2] = 0., b[2][0] = 0., b[2][1] = 0., b[2][2] = 0., c[1][0] = 0., c[1][1] = 0., c[1][2] = 0., c[2][0] = 0., c[2][1] = 0., c[2][2] = 0., e[1][0] = 0., e[1][1] = 0., e[1][2] = 0., e[2][0] = 0., e[2][1] = 0., e[2][2] = 0., g[1][0] = 0., g[1][1] = 0., g[1][2] = 0., g[2][0] = 0., g[2][1] = 0., g[2][2] = 0.}

[.5000000000 = .5, 0. = 0, 0. = 0, 0. = 0, 0. = 0, .5000000000 = .5, 0. = 0, 0. = 0, 0. = 0, 0. = 0, .5000000000 = .5, 0. = 0, 0. = 0, 0. = 0, 0. = 0, .5000000000 = .5, 0. = 0, 0. = 0, 0. = 0, 0. = 0, .5000000000 = .5, 0. = 0, 0. = 0, 0. = 0, 0. = 0, .5000000000 = .5, 0. = 0, 0. = 0, 0. = 0, 0. = 0]

@samen 

Error, (in fsolve) number of equations, 6, does not match number of variables, 18

 

 

Now do yo really believe that you can solve 6 equations  containing 18 unlnowns?? What planet are you from!???

So first of all, this error has to be fixed. You introduced semicolons between the equations within the for-loop. I mean - Why?

Now when I fix this basic problem, I get 30 equations in 30 unknowns - which is good!

But fsolve() still quits without a solution. This is a weakness in the fsolve() command: it is implementing an iterative process whilst trying to find a solution and "gives up" after a certain number of iterations. There is no user control over the number of iterations befoore it bails out.

Whay doess fsolve() quit? Well you have 30 equations in 30 unknowns, and if these were linear in the unlnowns, then absolutely no problem. However your system contains quations which are (at least) cubic in the unknowns. So 30 equations containing an unspecified amount of "cubics". Yep! that is going to get "difficult".

In the attached  have invoked a (free) Maple add-on package called DirectSearch() which can be somewhat more robust than the straightforward fsolve(), and certailnly allows more user control. This did find a "solution" once I adjusted one of its parameter options. The loctaed solution contains a maximum (absolute) residual error across all 30 equations of ~7e-04, and the sum of the squared residuals across all 30 equations is ~2e-06.

You will not be able to re-execute the attached until/unless you download/install the DirectSearch() package from the Maple Application centre here https://www.maplesoft.com/applications/view.aspx?SID=101333. Even once yu have this installed, you should be aware that it takes a few minutes to execute

restart; `ε` := 1.2; k := 2; M := 3

NULL

MyEqs := NULL; for i while i <= 2^(k-1)*M do t[i] := (i-.5)/(2^(k-1)*M); MyEqs := MyEqs, 9.797958972*a[1][1]+a[1][2]*(151.7893277*t[i]-37.94733192)+9.797958972*a[2][1]+a[2][2]*(151.7893277*t[i]-113.8419958)+(.3*(1.414213562*a[1][0]+a[1][1]*(9.797958972*t[i]-2.449489743)+a[1][2]*(4.743416490*(4*t[i]-1)^2-1.581138830)+1.414213562*a[2][0]+a[2][1]*(9.797958972*t[i]-7.348469229)+a[2][2]*(4.743416490*(4*t[i]-3)^2-1.581138830)))*(1.414213562*b[1][0]+b[1][1]*(9.797958972*t[i]-2.449489743)+b[1][2]*(4.743416490*(4*t[i]-1)^2-1.581138830)+1.414213562*b[2][0]+b[2][1]*(9.797958972*t[i]-7.348469229)+b[2][2]*(4.743416490*(4*t[i]-3)^2-1.581138830))-.2828427124*e[1][0]-.2*e[1][1]*(9.797958972*t[i]-2.449489743)-.2*e[1][2]*(4.743416490*(4*t[i]-1)^2-1.581138830)-.2828427124*e[2][0]-.2*e[2][1]*(9.797958972*t[i]-7.348469229)-.2*e[2][2]*(4.743416490*(4*t[i]-3)^2-1.581138830)+.9899494934*a[1][0]+.7*a[1][1]*(9.797958972*t[i]-2.449489743)+.7*a[1][2]*(4.743416490*(4*t[i]-1)^2-1.581138830)+.9899494934*a[2][0]+.7*a[2][1]*(9.797958972*t[i]-7.348469229)+.7*a[2][2]*(4.743416490*(4*t[i]-3)^2-1.581138830) = .5, 9.797958972*b[1][1]+b[1][2]*(151.7893277*t[i]-37.94733192)+9.797958972*b[2][1]+b[2][2]*(151.7893277*t[i]-113.8419958)-(.3*(1.414213562*a[1][0]+a[1][1]*(9.797958972*t[i]-2.449489743)+a[1][2]*(4.743416490*(4*t[i]-1)^2-1.581138830)+1.414213562*a[2][0]+a[2][1]*(9.797958972*t[i]-7.348469229)+a[2][2]*(4.743416490*(4*t[i]-3)^2-1.581138830)))*(1.414213562*b[1][0]+b[1][1]*(9.797958972*t[i]-2.449489743)+b[1][2]*(4.743416490*(4*t[i]-1)^2-1.581138830)+1.414213562*b[2][0]+b[2][1]*(9.797958972*t[i]-7.348469229)+b[2][2]*(4.743416490*(4*t[i]-3)^2-1.581138830))+(.3*(1.414213562*b[1][0]+b[1][1]*(9.797958972*t[i]-2.449489743)+b[1][2]*(4.743416490*(4*t[i]-1)^2-1.581138830)+1.414213562*b[2][0]+b[2][1]*(9.797958972*t[i]-7.348469229)+b[2][2]*(4.743416490*(4*t[i]-3)^2-1.581138830)))*(1.414213562*c[1][0]+c[1][1]*(9.797958972*t[i]-2.449489743)+c[1][2]*(4.743416490*(4*t[i]-1)^2-1.581138830)+1.414213562*c[2][0]+c[2][1]*(9.797958972*t[i]-7.348469229)+c[2][2]*(4.743416490*(4*t[i]-3)^2-1.581138830))+1.414213562*b[1][0]+1.0*b[1][1]*(9.797958972*t[i]-2.449489743)+1.0*b[1][2]*(4.743416490*(4*t[i]-1)^2-1.581138830)+1.414213562*b[2][0]+1.0*b[2][1]*(9.797958972*t[i]-7.348469229)+1.0*b[2][2]*(4.743416490*(4*t[i]-3)^2-1.581138830) = 0, 9.797958972*c[1][1]+c[1][2]*(151.7893277*t[i]-37.94733192)+9.797958972*c[2][1]+c[2][2]*(151.7893277*t[i]-113.8419958)-(.3*(1.414213562*b[1][0]+b[1][1]*(9.797958972*t[i]-2.449489743)+b[1][2]*(4.743416490*(4*t[i]-1)^2-1.581138830)+1.414213562*b[2][0]+b[2][1]*(9.797958972*t[i]-7.348469229)+b[2][2]*(4.743416490*(4*t[i]-3)^2-1.581138830)))*(1.414213562*c[1][0]+c[1][1]*(9.797958972*t[i]-2.449489743)+c[1][2]*(4.743416490*(4*t[i]-1)^2-1.581138830)+1.414213562*c[2][0]+c[2][1]*(9.797958972*t[i]-7.348469229)+c[2][2]*(4.743416490*(4*t[i]-3)^2-1.581138830))+3.676955261*c[1][0]+2.6*c[1][1]*(9.797958972*t[i]-2.449489743)+2.6*c[1][2]*(4.743416490*(4*t[i]-1)^2-1.581138830)+3.676955261*c[2][0]+2.6*c[2][1]*(9.797958972*t[i]-7.348469229)+2.6*c[2][2]*(4.743416490*(4*t[i]-3)^2-1.581138830) = 0, 9.797958972*e[1][1]+e[1][2]*(151.7893277*t[i]-37.94733192)+9.797958972*e[2][1]+e[2][2]*(151.7893277*t[i]-113.8419958)-.5656854248*c[1][0]-.4*c[1][1]*(9.797958972*t[i]-2.449489743)-.4*c[1][2]*(4.743416490*(4*t[i]-1)^2-1.581138830)-.5656854248*c[2][0]-.4*c[2][1]*(9.797958972*t[i]-7.348469229)-.4*c[2][2]*(4.743416490*(4*t[i]-3)^2-1.581138830)+2.969848480*e[1][0]+2.1*e[1][1]*(9.797958972*t[i]-2.449489743)+2.1*e[1][2]*(4.743416490*(4*t[i]-1)^2-1.581138830)+2.969848480*e[2][0]+2.1*e[2][1]*(9.797958972*t[i]-7.348469229)+2.1*e[2][2]*(4.743416490*(4*t[i]-3)^2-1.581138830) = 0, 9.797958972*g[1][1]+g[1][2]*(151.7893277*t[i]-37.94733192)+9.797958972*g[2][1]+g[2][2]*(151.7893277*t[i]-113.8419958)-1.697056274*e[1][0]-1.2*e[1][1]*(9.797958972*t[i]-2.449489743)-1.2*e[1][2]*(4.743416490*(4*t[i]-1)^2-1.581138830)-1.697056274*e[2][0]-1.2*e[2][1]*(9.797958972*t[i]-7.348469229)-1.2*e[2][2]*(4.743416490*(4*t[i]-3)^2-1.581138830)+.9899494934*g[1][0]+.7*g[1][1]*(9.797958972*t[i]-2.449489743)+.7*g[1][2]*(4.743416490*(4*t[i]-1)^2-1.581138830)+.9899494934*g[2][0]+.7*g[2][1]*(9.797958972*t[i]-7.348469229)+.7*g[2][2]*(4.743416490*(4*t[i]-3)^2-1.581138830) = 0 end do; numelems([MyEqs]); numelems(indets([MyEqs])); sol := DirectSearch:-SolveEquations([MyEqs], evaluationlimit = 1000000)

30

 

30

 

sol := [2.37592621476676*10^(-6), Vector(30, {(1) = HFloat(-2.342875231988728e-4), (2) = HFloat(6.890235854370985e-4), (3) = HFloat(9.555500400892925e-5), (4) = HFloat(-7.329665822908282e-5), (5) = HFloat(-8.162692029145546e-5), (6) = HFloat(-1.9859470194205642e-4), (7) = HFloat(-3.372940045665018e-5), (8) = HFloat(2.8733868020935915e-5), (9) = HFloat(2.3997301468625665e-5), (10) = HFloat(2.295104604854714e-4), (11) = HFloat(3.705027047544718e-4), (12) = HFloat(-7.278519369720016e-4), (13) = HFloat(-1.0384783490735572e-4), (14) = HFloat(8.33789526950568e-5), (15) = HFloat(-3.952924053010065e-4), (16) = HFloat(-8.985202293843031e-5), (17) = HFloat(-3.2907730565057136e-4), (18) = HFloat(-4.286115654394962e-6), (19) = HFloat(1.8699589418247342e-5), (20) = HFloat(3.876888040394988e-4), (21) = HFloat(4.47942328037243e-4), (22) = HFloat(-3.10515497403685e-5), (23) = HFloat(4.405116669659037e-5), (24) = HFloat(-9.557613429933554e-5), (25) = HFloat(-1.4027580255060457e-4), (26) = HFloat(-4.0726567385718226e-4), (27) = HFloat(4.7083772005862556e-4), (28) = HFloat(-7.78645699028857e-7), (29) = HFloat(4.733476089313626e-5), (30) = HFloat(1.5131048712646589e-5)}), [a[1][0] = -10058.4244568739, a[1][1] = -99527.1128306263, a[1][2] = 959.543739407914, a[2][0] = 341957.315489905, a[2][1] = 92094.5200906338, a[2][2] = -959.543750453665, b[1][0] = -22119.8252197932, b[1][1] = 38030.2869787544, b[1][2] = -6747.43057578954, b[2][0] = -19094.6538762761, b[2][1] = 14235.0833925461, b[2][2] = 6747.43071101906, c[1][0] = 4085.67720682402, c[1][1] = 6442.66289195766, c[1][2] = -2087.35427024835, c[2][0] = 1601.08016351940, c[2][1] = 9725.91375310733, c[2][2] = 2087.35427470547, e[1][0] = 9456.25823266217, e[1][1] = 4974.67327751617, e[1][2] = 1623.84535083997, e[2][0] = -48475.2035662693, e[2][1] = -17552.9252676696, e[2][2] = -1623.84535562594, g[1][0] = -45948.1353524541, g[1][1] = 69667.1727213737, g[1][2] = -12523.9743807409, g[2][0] = -27359.2840474511, g[2][1] = 27343.1159393932, g[2][2] = 12523.9743653372], 527688]

(1)

``

add(sol[2]^~2);max(sol[2]);min(sol[2]);

HFloat(2.3759262147667575e-6)

 

HFloat(6.890235854370985e-4)

 

HFloat(-7.278519369720016e-4)

(2)

``

Download DSAgain.mw

 

@Josolumoh 

The line

N=100

is an equation which happens to contain the symbol 'N'. It does not assign a value to 'N'. Using N at any point later in your worksheet (eg expecting it to evaluate to100) will fail miserably.

The line

X__1 := Statistics:-RandomVariable(Binomial(N, 4/10));

will initialise a random variable with thespecified statistic properties. No values of the random varibale will be generated, because you haven't asked for any. To get some values from this random variable, you will have to use the Statistics:-Sample() command.

If you really want to export to Excel, then the best way is to use the ExcelTools:-Export() command.

Fixing all of the above, the command sequence

  restart;
  N:=100;
  X__1 := Statistics:-RandomVariable(Binomial(N, 4/10));
  ExcelTools:-Export( Statistics:-Sample(X__1, 50),
                      "C:/Users/TomLeslie/Desktop/X1.xlsx"
                     );

will initialise a random variable with a distribution governed by Binomial(100, 4/10), and then export 50 samples of this random variable to the Excel file "C:/Users/TomLeslie/Desktop/X1.xlsx". Obviusly you will have to change this file path/name to something appropriate for your machine

@Gabriel samaila 

when you say

"even though is not specified at what value of phi(y) is given but i think is at phi(0). "

because in you original worksheet (and in all of the responses provided in this thread), you have specified the boundary conditions as

BCs := { phi(1) = 0, phi(eta) = 1, theta(1) = 0,
               theta(eta) = 1, v(1) = 0, v(eta) = 1
            }:

You have also hardwired eta=0.5, so soluions are only available for 0.5 <= y <= 1. Thus phi(0) is not defined by your ODE system and cannot be calculated

In the attached, I have

  1. Generated/plotted solutions of the ODE system for all of the values of 'Ha' and 'Sc' used in the table which you upplied as Sc.docx
  2. Generated an output table. for all of the values of 'Ha' and 'Sc' used in the table which you upplied as Sc.docx. Since you don't know what output function you want for this table, there is no way I can tell, so I have just used phi(y)*theta(y)/v(y) at y=0.75 as an illustrative example. You can get the value of any expression you want by reassigning the variable 'func' in the final execution group, and get it at any y-value you want (within the limits imposed by your boundary condition), by reassigning the variable 'argY'

  restart;
  kernelopts(version);

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

(1)

#
# Define ODEs
#
  ODEs:= {diff(phi(y), y, y)-Re*Sc*v(y)*(diff(phi(y), y))+Nt*(diff(theta(y), y, y)+(diff(theta(y), y))/y)/Nb = 0,
         diff(v(y), y, y)+(diff(v(y), y))/y-(Ha^2/(1-eta)^2+1/y^2)*v(y)-Re*v(y)*(diff(v(y), y)) = 0,
         diff(theta(y), y, y)+Ec*Pr*(diff(v(y), y)-v(y)/y)^2-Pr*Re*v(y)*(diff(theta(y), y))+Nb*(diff(theta(y), y))*(diff(phi(y), y))+Nt*(diff(theta(y), y))^2 = 0}:

#
# Define boundary conditions: note that these restrict
# the solution to the range eta < y < 1.0 (or possibly
# 1.0 < y < eta, depending on whether eta>1 or eta<1)
#
# In the parameter list provided by the OP, eta is
# hardwired to 0.5, so valid solution will be in the range
#
#  0.5 <= y <= 1.0
#
  BCs := { phi(1) = 0, phi(eta) = 1, theta(1) = 0,
           theta(eta) = 1, v(1) = 0, v(eta) = 1
         }:

#
# Define values for Ha and Sc which are used
# to generarte output table
#
  HaVals:=[0, 2, 4, 6]:
  ScVals:=[0.5, 1.0, 2.0]:
#
# Solve the ODE system for all pairs of values in
# the lists HaVals, ScVals, and park all solutions
# in a Maple Table()
#
  for j from 1 by 1 to numelems(HaVals) do
      for k from 1 by 1 to numelems(ScVals) do
        #
        # Generate parameter list
        #
          pVals:= [ eta = .5,       Ha = HaVals[j],
                    Sc = ScVals[k], Nt = .1,
                    Nb = .1,        Re = 2,
                    Ec = 0.1e-1,    Pr = 10
                  ];
        #
        # Produce ODE solution for this parameter list
        #
          sol[j,k]:=dsolve( `union`( eval(ODEs, pVals),
                                     eval(BCs, pVals)
                                   ),
                             numeric,
                             output=listprocedure
                           );
      od:
  od:

#
# generate plots for theta(y), phi(y), v(y) for
# all Ha and Sc values.
#
# Note that theta(y) and v(y) don't seem to vary
# as Sc changes value (but phi(y) does)
#
  pfunc:= z-> plot
              ( [ seq
                  ( seq
                    ( eval(z(y), sol[i,j])(y),
                      i=1..numelems(HaVals)
                    ),
                    j=1..numelems(ScVals)
                  ) ],
                  y=0.5..1,
                  title=cat( convert(z, string),
                             "(y) vs y"
                           ),
                  titlefont=[times, bold, 24],
                  legend=[ seq
                           ( seq
                             ( typeset
                               ( "Ha= ", HaVals[i],
                                 ", Sc= ", ScVals[j]
                               ),
                               i=1..numelems(HaVals)
                             ),
                             j=1..numelems(ScVals)
                           )
                         ],
                  legendstyle=[font=[times, bold, 16]],
                  size=[800, 600]
              ):
  pfunc(theta);
  pfunc(phi);
  pfunc(v);

 

 

 

#
# Define function for the output table. This could
# be any expression involving  phi(y), theta(y), v(y)
# so just for illustration, let's use
#
#   phi(y)*theta(y)/v(y)
#
  func:=phi(y)*theta(y)/v(y):
#
# Define the value of 'y' at which the output is
# required. The boundary conditions require that
# 0.5 < y < 1.0. For illustration, let's use
#
#   y=0.75
#
  argY:=0.75:
#
# Print header lines for output table
#
  printf
  ( "%30s\n %8s%12.6f%12.6f%12.6f\n\n",
    "Sc","Ha", seq
               ( ScVals[j],
                 j=1..numelems(ScVals)
               )
  );

  seq
  ( printf
    ( " %8d%12.6f%12.6f%12.6f\n%",
      HaVals[k],
      seq
      ( eval
        ( func, sol[k,j]
        )(argY),
        j=1..numelems(ScVals)
      )
    ),
    k=1..numelems(HaVals)
  );

                            Sc
       Ha    0.500000    1.000000    2.000000

        0    0.165945    0.248300    0.408680
        2    0.383798    0.462637    0.616654
        4    1.396358    1.481539    1.647812
        6    4.323298    4.439670    4.666815

 

 

Download odeProb6.mw

First 72 73 74 75 76 77 78 Last Page 74 of 207