Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 163 days
Contact tomleslie

MaplePrimes Activity


These are answers submitted by tomleslie

Crude way...

tomleslie 13876
October 07 2021
2 0

try using 100 raised to the power o (as in the letter 'o') ie 100^o

Obvious question...

tomleslie 13876
October 06 2021
1 0

What do you want to happen when more than one element satisfied the criterion??

While you think about that, you might consider the "toy" example shown in the attached whch shows a couple of different ways

> 

  restart;
  L:=[seq( j/2, j=10..0,-1)];
  seq( `if`( `and`(L[j]>0, L[j]<1), [j,L[j]], NULL), j=1..numelems(L));
  ListTools:-SelectFirst(numelems(L), j-> `and`(j>0, j<1), L, output=[indices, values]);

[5, 9/2, 4, 7/2, 3, 5/2, 2, 3/2, 1, 1/2, 0]

  

[10, 1/2]

  

[10], [1/2]

 (1)
> 

 

Download listSel.mw

In Maple speak...

tomleslie 13876
October 05 2021
1 1

This is known as the "Context Bar"

So make sure it is toggled in the View menu

One way...

tomleslie 13876
October 05 2021
1 0

is shown in the attached. (The plot renders rather differently in the worksheet thsn it does on this site!)

> 

  restart;
  with(Statistics):
  x:= RandomVariable(Binomial(45, 0.9)):
  data:=[seq(ProbabilityFunction(x, i), i=0..45)]:
  dataplot(data);

 
> 

 

Download plFata.mw

Answers are the same...

tomleslie 13876
October 01 2021
2 4

When you compute the "same" quantity in two completely different ways, you really shouldn't be surprised if the answers "look" different. The important question is - are they different?

The attached shows that both methods  in fact give the same answer, after appropriate simplification. Check the final entry from each execution group.

> 

  restart:
  with(Syrup):
  with(inttrans):
#
# Set up the problem
#
  ckt:= [ vin(5), Rm(1.00000000*10^7), Cm(1.80000000*10^(-11), ic = 0),
          R1(9.00000000*10^6) &// C1(1.50000000*10^(-11), ic = 0),
          R2(1.00000000*10^6) &// C2(1.50000000*10^(-11), ic = 0)
        ]:
  Draw(ckt);

 
> 

#
# Get the voltage across the capacitor by solving
# the odes (basically a transient analysis)
#
  vcm:= eval
        ( v[Cm](t),
          dsolve
          ( Solve
            ( ckt,
              'tran'
            )
          )
        );
  simplify(expand(vcm));

  

(1/3)*exp((100000/459)*(-119+8041^(1/2))*t)*((61/64328)*8041^(1/2)-1/8)*8041^(1/2)-(1/3)*exp(-(100000/459)*(119+8041^(1/2))*t)*(-(61/64328)*8041^(1/2)-1/8)*8041^(1/2)+(91/3)*exp((100000/459)*(-119+8041^(1/2))*t)*((61/64328)*8041^(1/2)-1/8)+(91/3)*exp(-(100000/459)*(119+8041^(1/2))*t)*(-(61/64328)*8041^(1/2)-1/8)+5/2

  

(1/32164)*(-415*8041^(1/2)-40205)*exp((100000/459)*(-119+8041^(1/2))*t)+5/2+(1/32164)*(415*8041^(1/2)-40205)*exp(-(100000/459)*(119+8041^(1/2))*t)

 (1)
> 

#
# Get the transfer function to the node '2' (basically
# an ac analysis) then compute the step response, and
# convert the result from hyperbolic trig functions to
# exponentials
#
  v2:= convert
       ( invlaplace
         ( eval
           ( v[2],
             [ Solve
               ( ckt,
                 returnall
               )
             ][1]
           )/s,
           s,
           t
        ),
        exp
      );
  simplify(expand(v2));

5/2-(5/32164)*(83*exp((200000/459)*t*8041^(1/2))*8041^(1/2)+8041*exp((200000/459)*t*8041^(1/2))-83*8041^(1/2)+8041)*exp(-(700000/27)*t-(100000/459)*t*8041^(1/2))

  

(1/32164)*(-415*8041^(1/2)-40205)*exp((100000/459)*(-119+8041^(1/2))*t)+5/2+(1/32164)*(415*8041^(1/2)-40205)*exp(-(100000/459)*(119+8041^(1/2))*t)

 (2)
> 

 

Download ckt.mw

|Several ways to do this...

tomleslie 13876
September 30 2021
3 5

of which one is shown in the attached

> 

  restart;
  with(RootFinding):
  sols:= []:
  tmp:=0:
  while true do
        tmp:= NextZero
              (  x -> sin(x)-cos(4*x - 1/6*Pi),
                 tmp
              );
        if   tmp <= evalf(2*Pi) then
             sols:= [sols[], tmp];
        else break;
        fi:
  end do:
  sols;

[.4188790204, 1.675516081, 1.745329251, 2.932153143, 3.839724354, 4.188790204, 5.445427266, 5.934119456]

 (1)
> 

 

Download roots.mw

 

Set ordering...

tomleslie 13876
September 27 2021
1 1

You may have to use Google Translate!

There is no guaranteed order in a set. Mathematically, the position of a particular element in a set has no meaning.

As a practical matter, Maple will generally group set elements by their type, and then alphanumerically within each type.

Consider the "toy" example attached where a set (where order does not matter) is constructed from a list (where order does matter). Note how the "order" changes

> 

L:=[y(x), z(x),3,2,1,b,c,a];
S:={}:
for i from 1 by 1 to numelems(L) do
    S:= S union {L[i]};
od:
S;

[y(x), z(x), 3, 2, 1, b, c, a]

  

{1, 2, 3, a, b, c, y(x), z(x)}

 (1)
> 

 

Download setList.mw

Missing multiplication signs...

tomleslie 13876
September 27 2021
1 0

See the annotated corrections in the attached

> 

 omega1:=c^2/2+1/2*sqrt(Ds);
 omega2:=c^2/2-1/2*sqrt(Ds);

(1/2)*c^2+(1/2)*Ds^(1/2)

  

(1/2)*c^2-(1/2)*Ds^(1/2)

 (1)
> 

 Ds:=(Ea^2-Eb^2)^2+4*Eab*X5^2;

(Ea^2-Eb^2)^2+4*Eab*X5^2

 (2)
> 

 Ea:=sqrt((epsk1+X1)*(epsk2+X2)); # <-missing multiplcation sign
 Eb:=sqrt((epsk2+X2)*(epsk2+X1)); # <-missing multiplcation sign
 Eab:=sqrt((epsk1+X2)*(epsk2+X1)); # <-missing multiplcation sign

((epsk1+X1)*(epsk2+X2))^(1/2)

  

((epsk2+X2)*(epsk2+X1))^(1/2)

  

((epsk1+X2)*(epsk2+X1))^(1/2)

 (3)
> 

 c:=sqrt(Ea^2+Eb^2);
 epsk1:=k^2/2/m1;
 epsk2:=k^2/2/m2;

((epsk1+X1)*(epsk2+X2)+(epsk2+X2)*(epsk2+X1))^(1/2)

  

(1/2)*k^2/m1

  

(1/2)*k^2/m2

 (4)
> 

 Ek4:=sqrt((Ea^2*Eb^2-Eab^2*X5^2)*Ds);

((((1/2)*k^2/m1+X1)*((1/2)*k^2/m2+X2)^2*((1/2)*k^2/m2+X1)-((1/2)*k^2/m1+X2)*((1/2)*k^2/m2+X1)*X5^2)*((((1/2)*k^2/m1+X1)*((1/2)*k^2/m2+X2)-((1/2)*k^2/m2+X2)*((1/2)*k^2/m2+X1))^2+4*(((1/2)*k^2/m1+X2)*((1/2)*k^2/m2+X1))^(1/2)*X5^2))^(1/2)

 (5)
> 

series(Ek4,k=0,4);

series(2*(X1*X2*(X1*X2-X5^2)*(X1*X2)^(1/2)*X5^2)^(1/2)+((1/4)*(X1*X2*(X1*X2-X5^2)*(X1*X2)^(1/2)*X5^2)^(1/2)*(4*X1^2*X2*m1+X1^2*X2*m2+3*X1*X2^2*m1+2*X1*X2^2*m2-3*X1*X5^2*m2-3*X2*X5^2*m1)/(m1*m2*X1*X2*(X1*X2-X5^2)))*k^2+O(k^4),k,4)

 (6)
> 

series(sqrt(Ek4),k=0,4);

series(4^(1/4)*(X1*X2*(X1*X2-X5^2)*(X1*X2)^(1/2)*X5^2)^(1/4)+((1/16)*4^(1/4)*(X1*X2*(X1*X2-X5^2)*(X1*X2)^(1/2)*X5^2)^(1/4)*(4*X1^2*X2*m1+X1^2*X2*m2+3*X1*X2^2*m1+2*X1*X2^2*m2-3*X1*X5^2*m2-3*X2*X5^2*m1)/(m1*m2*X1*X2*(X1*X2-X5^2)))*k^2+O(k^4),k,4)

 (7)
> 

 

Download missMult.mw

Is the function call necessary?...

tomleslie 13876
September 26 2021
1 0

Seems like overkill?. See the attached

> 

   restart;
#
# Uncomment the randomize() command below only if you
# want to generate a different set of "random" numbers
# each time the worksheet is executed - otherwise the
# same sequence of random numbers will be generated -
# so maybe not very(?) random
#
# Uncommenting the randomize() command is useful if
# you want to be *really* random, but it can make debug
# of subsequent operations difficult because, if an error
# occurs anywhere, you may not be able to reproduce it
#
#
# randomize():
#
   A := [1,1,2,5,6]:
   r := rand(1..numelems(A)):
   A[r()];
   A[r()];
   A[r()];
   A[r()];
   A[r()];
   A[r()];
   A[r()];
   A[r()];
   A[r()];
   A[r()];
   A[r()];
   A[r()];
   A[r()];
   A[r()];
   A[r()];
   A[r()];

6

  

1

  

6

  

1

  

2

  

5

  

5

  

6

  

2

  

1

  

6

  

1

  

2

  

1

  

1

  

5

 (1)
> 

 

Download rand.mw

 

Maybe an animation...

tomleslie 13876
September 26 2021
1 3

will suffice - see the attached

``

> 

P1 := plot(5*sin(3*theta), theta = 0 .. 2*Pi, coords = polar, color = magenta, thickness = 3, filled = [color = green, transparency = .5])

 
> 

P2 := plot(5*cos(3*theta), theta = 0 .. 2*Pi, coords = polar, color = red, thickness = 3, filled = [color = blue, transparency = .5])

 
> 

P3 := plot(5, theta = 0 .. 2*Pi, coords = polar, color = black, thickness = 3, filled = [color = yellow, transparency = .9])

 
> 

plots:-display(P1, P2, P3, axes = none)

 
> 

plots:-animate(plot, [[5*sin(a*theta), 5*cos(a*theta), 5], theta = 0 .. 2*Pi, coords = polar, color = [green, red, black], thickness = 3, filled = [color = blue, transparency = .5]], a = 2 .. 6)

 

NULL

NULL


 

Download anim2.mw

Some observations...

tomleslie 13876
September 26 2021
1 4

I agree with Preben that the function PhiB() is undefined. The function Ub() also appears to be undefined.

ode1 is a first-order equation - so can only support one initial (boundary??!!) condition. Most people would consider this to be an IVP rather than a BVP.

It seems superfluous to apply a "shooting" method to such a system. It is technically possible - but why?

It is possible to solve the equation ode1, subject to one initial condition: I chose phi(0)=0.

Since the original ODE is of third degree, multiple solutions are possible depending on which roots are selected. In fact nine solutions occur. However eight of these produce complex values which I assume(?) you don't want

The attached shows a plot of the only real solution

> 

  restart:
  dpdexi := 2:
  K := 0.2:
  epsilon := 0.0001:
  ode1 := (diff(phi(eta), eta) + epsilon*K*diff(phi(eta), eta)^3)/(1 + epsilon*diff(phi(eta), eta)^2) = dpdexi*eta:
  bc1 := phi(0) = 0, phi(1) = 0:
 
> 

#
# Solve the ode with a single initial condition
#
  dsolve([ode1, phi(0)=0]):
#
# Get all the possible values from RootOfs
#
  sol:=[allvalues(%)]:
#
# How many different solutions are there?
#
  numelems(sol);
#
# Hmmm 9 solutions, But only one of these is real
#
  solR:=[seq( `if`( type(simplify(rhs(evalf(eval(j, eta=1/2))), zero), realcons),j, NULL), j in sol)]:
  numelems(solR);
#
# Plot the real solution
#
  plot( rhs(solR[]), eta=0..1);

9

  

1

  

 
> 

 

Download odeProb.mw

 

Well.......

tomleslie 13876
September 25 2021
2 0

if you just want to replicate the plot - then you are certaily doing it the hard way!

See the attached which just produce the plot from the function/parameters in the image you provided - ie it makes no attempt to "derive" the relevant function.

> 

  restart;
  with(plots):
  params:= [ sigma=4,
             h=1,
             M=1,
             P=1
           ]:
  local Psi:
  Psi:= unapply
        ( eval
          ( exp( (-1/2)*x^2/(sigma^2+I*h*t/M)+(I*P*x/h-I*P^2/(2*M)*(t/h))/(1+I*h*t/(M+sigma^2))),
            params
          ),
          [x,t]
        );
  animate
  ( plot,
    [ [ abs(Psi(x,t)),
        Re(Psi(x,t)),
        Im(Psi(x,t))
      ],
      x=-100..100,
      color=[black, blue, red]
    ],
    t=0..40
  );

proc (x, t) options operator, arrow; exp(-(1/2)*x^2/(16+I*t)+(I*x-((1/2)*I)*t)/(1+((1/17)*I)*t)) end proc

  

 
> 

 

> 

 

> 

 


 

Download anim.mw

Hmmm??...

tomleslie 13876
September 24 2021
1 0

This code appears to work??????

So what exactly is the problem?

See the attached

> 

restart
interface(version)

`Standard Worksheet Interface, Maple 2021.1, Windows 7, May 19 2021 Build ID 1539851`

 (1)
> 

f := proc (t, x) options operator, arrow; 1+t*sin(x) end proc; t[0] := 0; x[0] := 0; h := .1

proc (t, x) options operator, arrow; 1+t*sin(x) end proc

  

0

  

0

  

.1

 (2)
> 

for n to 20 do t[n] := n*h; m1 := f(t[n-1], x[n-1]); m2 := f(t[n-1]+(1/2)*h, x[n-1]+(1/2)*m1*h); m3 := f(t[n-1]+(1/2)*h, x[n-1]+(1/2)*m2*h); m4 := f(t[n-1]+h, h*m3+x[n-1]); x[n] := x[n-1]+(1/6)*h*(m1+2*m2+2*m3+m4) end do
> 

data1 := [seq([t[n], x[n]], n = 0 .. 20)]; plot(data1, t = 0 .. 2, x = 0 .. 19)

 
> 

NULL
> 

NULL


 

Download grpplt.mw

 

Use the parameterranges option...

tomleslie 13876
September 22 2021
1 0

as in the "toy" example below

  with(Statistics):
  with(plots):
  X := Vector([1, 2, 3, 4, 5, 6], datatype=float):
  Y := Vector([2.2, 3, 4.8, 10.2, 24.5, 75.0], datatype=float):
  fit1:=NonlinearFit(a+b*v+exp(c*v), X, Y, v, parameterranges=[a=6..7, b=-5..-4, c=0..1]);
  fit2:=NonlinearFit(a+b*v+exp(c*v), X, Y, v, parameterranges=[a=7..8, b=-5..-4, c=0..1]);
  p1:=plot(X,Y, color=red, style=point, symbol=solidcircle, symbolsize=20):
  p2:=plot(fit1, v=1..6, color=green):
  p3:=plot(fit2, v=1..6, color=blue):
  display([p1, p2, p3]);

fit1 in the above is an "accurate" fit. Note that in fit2, one of the parameters is right at the end of its allowed range

Recursive definition...

tomleslie 13876
September 20 2021
1 0

You have the defined procedure

U := proc(k, h)
                options operator, arrow;
                (sum(sum(U(r, h-s)*U(k-r, s), s = 0 .. h), r = 0 .. k)+(k+1)*U(k, h))/m
           end proc;

Note that this procedure is assigned to the name 'U', but contains calls to a procedures also called U() - in other words the procedure calls itself. Such a recursive definition can be useful in certain programs, but always needs  a stopping criterion - otherwise the recursion goes on for ever!

First 40 41 42 43 44 45 46 Last Page 42 of 207