mmcdara

3437 Reputation

16 Badges

5 years, 230 days

MaplePrimes Activity


These are answers submitted by mmcdara

Did you notice that sec(theta - (Pi/2) * (2/Pi *(theta + Pi/4))) = sqrt(2) ?

sec.mw


 

restart:
eq := x^6+x^4*(1.02-0.50e-2*y)+x^2*(-0.4950e-2*y+0.625e-5*y^2)+0.625e-5*y^2=0:

s := [solve(eq, x)]:
# check the number of solutions
numelems(s);
                               6

# define a function (sol) of y onto s
sol := unapply(s, y):

# Arrange the 6 plots of real and imaginay parts of each solution 
# into a 2x3 matrix 

M := Matrix(
      2, 3,
      (i, j) -> plot([(Re, Im)(sol(y)[3*(i-1)+j])], y=-10..10, color=[blue, red], legend=['Re', 'Im'], title=cat("Solution ", 3*(i-1)+j))
     ):

# Display M
DocumentTools:-Tabulate(M, width=60);

A_solution.mw

You will find in the first part (2D input mode) of the attached file something which runs without error.

But the results are very strange and I suspect something's wrong in your FDTD scheme.

So the last part of this same file (1D input mode) contains a much simpler way to obtain the 5 points discretization of the laplacian (matrix form) and a time marching based on a simple forward euler scheme (dt=0.1).
The Explore command (could have been animate instead) will help you to verify qualitatively that the diffusion process seemsto perform correctly.
The peak at x=y=0 is the result of he overlap of the BCs.at this point (given the additive way you constructed your BCs you have counted twice the boundary corner (0,0)).

Download 2d_heat_eq_RK4_mmcdara.mw

Adaptation to a RK4 time scheme is easy

If you just want the extrema ( @vv I didn't know about the term exprema) you can also use this
(this method is slower than the classical vv's method but, at least in this case, it works well and requires almost no mathematical background)

restart:
with(Student[Calculus1]):
A := x^2 + 400/x:
ExtremePoints( A );
                           [   (2/3)]
                           [2 5     ]

See also the attached file
AAA.mw


In fact the true question is "What form would you like this sensitivity report to take?"

Here are a few solutions (I used a larger step to provide "small" displays.
The result matrix can be saved using ExportMatrix for instance.

SR.mw

Unless I made a mistake here is a proof
 

restart:

g := exp(x)+exp(y)+exp(z)+2*exp(-x-y-z)-4*2^(1/4);

exp(x)+exp(y)+exp(z)+2*exp(-x-y-z)-4*2^(1/4)

(1)

# Write
#   u=exp(x)
#   v=exp(y)
#   w=exp(z)

h := u+v+w+2/(u*v*w)-4*2^(1/4)

u+v+w+2/(u*v*w)-4*2^(1/4)

(2)

# Observe that:
#   u+v+w is 3 times the Arithmetic Mean (AM) of (u, v, w)
#   u*v*w equals the Geometric Mean (GM) of (u, v, w) raised to the cube
# Thus

h := 3*AM+2/GM^3-4*2^(1/4)

3*AM+2/GM^3-4*2^(1/4)

(3)

# Let's focus on the term 3*AM+2/GM^3
#
# From the AM–GM inequality
# https://en.wikipedia.org/wiki/Inequality_of_arithmetic_and_geometric_means
# one has GM <= AM
# Thus:
#   GM^3 <= AM^3
#   1/GM^3 >= 1/AM^3
#   3*AM+2/GM^3 >= 3*AM+2/AM^3
#
# Solve the equation 3*AM+2/AM^3 > 4*2^(1/4) given that AM > 0

solve(3*AM+2/AM^3 > 4*2^(1/4)) assuming AM > 0

RealRange(Open(0), Open(2^(1/4))), RealRange(Open(2^(1/4)), infinity)

(4)

# Thus 3*AM+2/AM^3 is
#    positive provided AM <>2^1/4,
#    null if AM = 2^1/4,
# and, because AM is necessarilly > 0 (arithemic mean of exponentials),
# this proves 3*AM+2/AM^3 > 4*2^(1/4) is true but at AM = 2^1/4.
#
# Given that 3*AM+2/GM^3 >= 3*AM+2/AM^3 you have your proof.

 

Download AM_GM.mw

An even simpler way based on th same AM-GM inequality is
 

# Other way: AM = GM iif u=v=w.
# Then it's enough to demonstrate that 3*u+2/u^3 cannot be less than 4*2^(1/4)
# for any strictly positive u:

infolevel[solve]:=1:
solve({3*u+2/u^3 < 4*2^(1/4), u > 0});
solve: Warning: no solutions found

 

Possibly apart the obvious loops at edges x0 and x10(*) I can't see any cycle (more precisely circuit as your graph is directed) on your graph !
A cycle circuit in a directed graph G is a non-empty path in G with its first vertex equal to its last vertex.

(*) some will consider that a loop is not circuit; ths choice to do this or not generally depends upon the context

If I correctly understood your problem you will find several ways to do this in the attached file.
 

restart:

One way is to define a function with 4 indeterminates

f   := (a, b, c, x )-> a*x^2+b*x+c;
sol := solve(f(a, b, c, x)=0, x);

# ... the following is strictly identical to what id done below.

 

proc (a, b, c, x) options operator, arrow; a*x^2+b*x+c end proc

 

(1/2)*(-b+(-4*a*c+b^2)^(1/2))/a, -(1/2)*(b+(-4*a*c+b^2)^(1/2))/a

(1)

Another way is to define a parameterized function of x this way

f := (a, b, c) -> x -> a*x^2+b*x+c;

proc (a, b, c) options operator, arrow; proc (x) options operator, arrow; a*x^2+b*x+c end proc end proc

(2)

sol   := solve(f(a, b, c)(x)=0, x):
optx  := unapply( sol[1], (a, b, c));
opttx := unapply( sol[2], (a, b, c));
 

proc (a, b, c) options operator, arrow; (1/2)*(-b+(-4*a*c+b^2)^(1/2))/a end proc

 

proc (a, b, c) options operator, arrow; -(1/2)*(b+(-4*a*c+b^2)^(1/2))/a end proc

(3)

# to evaluate optx at b=2 and c=4 just do

optx(a, 2, 3);
 

(1/2)*(-2+(-12*a+4)^(1/2))/a

(4)

# be careful, this value can be complex

eval(optx(a, 2, 3), a=-1);
evalf(%);

eval(optx(a, 2, 3), a=1);
evalf(%);

1-(1/2)*16^(1/2)

 

-1.000000000

 

-1+(1/2)*(-8)^(1/2)

 

-1.+1.414213562*I

(5)

# A plot without caution (Maple 2015)
# Note that the plot is right for a=-1, but wrong for a=+1

plot(optx(a, 2, 3), a=-2..2, gridlines=true);

 

# A more convenient way

plot([Re(optx(a, 2, 3)), Im(optx(a, 2, 3))], a=-2..2, color=[blue, red], legend=['Re(optx)', 'Im(optx)'], gridlines=true);

 

A third way could be to define this parameterized function as below, but it seems less natural to me
(my proper opinion)

f   := x -> (a, b, c) -> a*x^2+b*x+c;
sol := solve(f(x)(a, b, c)=0, x);

proc (x) options operator, arrow; proc (a, b, c) options operator, arrow; a*x^2+b*x+c end proc end proc

 

(1/2)*(-b+(-4*a*c+b^2)^(1/2))/a, -(1/2)*(b+(-4*a*c+b^2)^(1/2))/a

(6)

 


 

Download Several_ways.mw

This just says that the solution of F(T, A)=0 with respect to T expresses itself in terms of the LambertW function.

Your solution cannot be expressed in another way given the definition of the LambertW function.
(In fact this latter has been introduced to designate the series solution of the the equation you want to solve; see
here Lambert_W_function (History) for more details)
 

restart

df := (-7.233333422*10^10*exp(.45*T)+1.617777797*10^11)/T-(-1.607407427*10^11*exp(.45*T)+1.617777797*10^11*T+1.745679033*10^11+A)/T^2

(-0.7233333422e11*exp(.45*T)+0.1617777797e12)/T-(-0.1607407427e12*exp(.45*T)+0.1617777797e12*T+0.1745679033e12+A)/T^2

(1)

solve(df, T);

2.222222222*LambertW(-.3995249844-0.2288650873e-11*A)+2.222222222

(2)

# What does LambertW mean?

#help(LambertW)

# Example

exple := (a*x+b)*exp(k*x)=c;;
solve(exple, x);

(a*x+b)*exp(k*x) = c

 

(a*LambertW(c*k*exp(b*k/a)/a)-b*k)/(a*k)

(3)

# Assuming T <> 0, df=0 if numer(df)=0

ndf := numer(df);

-0.7233333422e11*T*exp(.45*T)-0.1745679033e12+0.1607407427e12*exp(.45*T)-A

(4)

# Put this into a more convenient form

expr := subs(exp(.45*T)=freeze(exp(.45*T)), ndf):
e    := indets(%, name)[-1]:
conv := thaw(collect(expr, e))

(-0.7233333422e11*T+0.1607407427e12)*exp(.45*T)-0.1745679033e12-1.*A

(5)

# "conv" has the same form of "exple" with
#  a = -7.233333422*10^10
#  b = 1.607407427*10^11
#  c = 1.745679033*10^11+1.*A
#  k = .45
#
# thus the solution you get.

 


 

Download LambertW.mw

 


The term probability function for a discrete random variable (RV)  should be conveniently replaced by the term probability mass function (PMF).
The PMF of a discrete RV defined over a countable set S of integers (for instance) is null evertwhere but at the points in S.

I consider IMO that the plots below are more correct.

 

restart;

with(Statistics):

X := RandomVariable(Binomial(45, 0.9)):
P := seq(ProbabilityFunction(X, i), i=0..45):

plots:-display(
  seq(
    plot([ [i, 0], [i, P[i+1]] ], thickness=5),
    i=1..45
  )
);

 

reduced_P := map(n -> if P[n+1] > max([P])/1000 then
                    plot([ [n, 0], [n, P[n+1]] ], thickness=7)
                    end if,
                    [$0..45]
                ):
plots:-display(reduced_P)

 

 


Download Probability_plot.mw


PS,
Let X a discrete RV.
Still IMO I don't agree with Maple returning a non zero value for ProbabilityFunction(X, t) when t doesn't belong to set S.
If this results was good, then the CDF of X wouldn't be a constant in the interval (s[i], s[i+1]) where s[i] and s[i+1] are two consecutive points in S (nor even a discontinuous function).
Use this to convince yourself that Maple is not consistent with the definitions it uses.

# Integrate the probability function between 0 and s
# One might expect this is the CDF of X

F := int(ProbabilityFunction(X, t), t=0..s);

# Compare the plots of F and ths CDF of X

plots:-display(
  plot(F, s=30..46, color=red),
  plot(CDF(X, s), s=30..46, color=blue)
)

 

 This generates a random polynomial whose coefficients are randomly choosen in the list L.

L := [1, 1, 2, 5, 6]:
randpoly(z, coeffs = proc() combinat:-randperm(L)[1] end proc);

 


EDITED 

I hadn't correctly understood your question, sorry.
There is no difficilty at all to do this with Explore:

 

f := proc(a)
plots:-display(
  plot(5, theta = 0 .. 2*Pi, coords = polar, color = black, thickness = 3, filled = [color = yellow, transparency = .9]),
  plot(
    5*sin(a*theta), 
    theta = 0 .. 2*Pi, 
    coords = polar, 
    color = magenta, 
    thickness = 3, 
    filled = [color = green, transparency = .5]
  ),
  plot(
    5*cos(a*theta), 
    theta = 0 .. 2*Pi, 
    coords = polar, 
    color = red, 
    thickness = 3, 
    filled = [color = blue, transparency = .5]
  )
)
end proc:


Explore(f(a), parameters = [a = 2 .. 6])


Download PétalesGraphes_Maple1015.mw

Point 1
What quantity has dimension 1/time^2 : no one, but it shoulf not be a problem
(of course an acceleration divided by a mass is homogenous tio the inverse of  time squared but I don't think its the answer you wnanted).

Point 2
By integration both sides of DIV from 0 and 3 one can avoid computing an integral to find the length of the wire:

a:=0.36: 
DIV := diff(y(x), x, x) = a*sqrt(1 + diff(y(x), x)^2):
RV  := y(0)=2.6, y(3)=2.1:
sol := rhs(dsolve({DIV, RV}, y(x))):
L   := evalf( ( eval(diff(sol, x), x=3) - eval(diff(sol, x), x=0) ) / a):
evalf(value(%));

                          3.187401749


 

restart:

DIV := diff(y(x), x, x) = a*sqrt(1 + diff(y(x), x)^2):
RV  := y(0)=2.6, y(3)=2.1:
sol := rhs(dsolve({DIV, RV}, y(x))):

L    := evalf( ( eval(diff(sol, x), x=3) - eval(diff(sol, x), x=0) ) / a):
vals := evalf([allvalues(L)]):  # 2 opposite solutions
sag  := unapply(vals, a):

plot(
  [abs(sag(a)[1]), abs(sag(a)[2])],
  a=0.2..0.8,
  color=[red, blue],
  thickness=[1, 11],
  transparency=[0, 0.8],
  gridlines=true,
  labels=['a', 'length']
)

 

# which value of a gives a length of 3.6?

rel := a = select(is, fsolve~(vals=~3.6), positive)[]

a = .6900379205

(1)

# Value of the Sag defined as the maximum vertical distance
# between the shape of the wire and the straitght line
# which passes by through the points (0, 2.6) and (3, 3.1):

sola  := eval(sol, rel):
vdist := sola - (2.6+(2.1-2.6)*x/3):

Sag := -minimize(vdist, x=0..3, location);

.8569225658, -{[{x = 0.}, -0.80202e-8], [{x = 1.462158429}, -.8569225658], [{x = 3.}, -0.8651e-8]}

(2)

# thus the solution

Sag__xy := op(2, [Sag][2])[2]

[{x = 1.462158429}, -.8569225658]

(3)

 


Download catenary_4.mw

Alternative to minimize:

x = select(is, [solve(diff(vdist, x)=0, x)], positive)[]:
allvalues([eval(vdist, %)]):
y = select(is, op~([%]), real)[]:

Sag__xy= [%%%, %];
         Sag__xy = [x = 1.462158429, y = -0.8569225635]

 

As Ì said in my answer to tyour previous question, this can be done in a simpler way

restart;
# if you use integer values, I advice you to scale the equations 
# by using a scaling factor C (this will help fsolve finding the solution).
# If you use floats instead this is not necessary

C := 100:

DIV := diff(y(x), x, x) = a*sqrt(1 + diff(y(x), x)^2):
RV  := y(0)=50/C, y(d)=70/C:

sol := rhs(dsolve({DIV, RV}, y(x))):
L   := Int(sqrt(1+diff(sol, x)^2), x=0..d) assuming x>0, d>0:

# abscissa of the contact to the ground

c := solve(diff(sol, x)=0, x):
# equations to find a and d:

eq1 := eval(sol, x=c)=0:
eq2 := value(L)=140/C:
ad  := fsolve({eq1, eq2}, {a, d})

              {a = 9.213375344, d = 0.5541723235}

# unscaled a and d

AD := [a, d]=~eval([a/C, d*C], ad)
              [a = 0.09213375344, d = 55.41723235]

 

@brian bovril 

The boundary conditions you used are not correct; moreover @vv  solution wich consists in solving two half wires seems a little bit artificial and complicated (IMO) .
If the wire has no transverse stiffness (whch is the basic assumptions used to derive the catenary equation), 
then the boundary conditions are 

y(0) = h0, y(d)=h1

where d is the distance between the two poles.

Here is a partial numeric solution for two poles of heights 1 and 2 and a cable of length 4 (just an example to adapt to the true values [see foot note])
 

restart;

# A simpler problem:
#
# Find the distance between two poles of height 1 and 2 such
# that the wire has length 4.

DIV := diff(y(x), x, x) = sqrt(1 + diff(y(x), x)^2):
RV  := y(0)=1, y(d)=2:
sol := rhs(dsolve({DIV, RV}, y(x))):
L   := Int(sqrt(1+diff(sol, x)^2), x=0..d) assuming x>0, d>0:
gap := solve(value(L)=4, d);

plot(eval(sol, d=gap), x=0..gap)

ln(17/2+(1/2)*285^(1/2))

 

 

# As this problem has a single solution we need an extra parameter
# to find this more complicated problem:
#
# Find the distance between two poles of height 1 and 2 such
# that the wire has length 4 AND the wire tangents the ground.

DIV := diff(y(x), x, x) = a*sqrt(1 + diff(y(x), x)^2):
RV  := y(0)=1, y(d)=2:
sol := rhs(dsolve({DIV, RV}, y(x)));
L   := Int(sqrt(1+diff(sol, x)^2), x=0..d) assuming x>0, d>0;

cosh(a*x+ln(RootOf(((exp(a*d))^2-exp(a*d))*_Z^2-2*a*exp(a*d)*_Z-exp(a*d)+1)))/a+(1/2)*(2*a*(exp(a*d))^2*RootOf(((exp(a*d))^2-exp(a*d))*_Z^2-2*a*exp(a*d)*_Z-exp(a*d)+1)-4*a*exp(a*d)*RootOf(((exp(a*d))^2-exp(a*d))*_Z^2-2*a*exp(a*d)*_Z-exp(a*d)+1)-(exp(a*d))^2+1)/(exp(a*d)*(exp(a*d)-1)*a*RootOf(((exp(a*d))^2-exp(a*d))*_Z^2-2*a*exp(a*d)*_Z-exp(a*d)+1))

 

Int((1+sinh(a*x+ln(RootOf(((exp(a*d))^2-exp(a*d))*_Z^2-2*a*exp(a*d)*_Z-exp(a*d)+1)))^2)^(1/2), x = 0 .. d)

(1)

# abscissa of the contact to the ground

c := solve(diff(sol, x)=0, x)

-ln(RootOf(((exp(a*d))^2-exp(a*d))*_Z^2-2*a*exp(a*d)*_Z-exp(a*d)+1))/a

(2)

# equations to find a and d:

eq1 := eval(sol, x=c)=0:
eq2 := value(L)=4:
ad  := fsolve({eq1, eq2}, {a, d})

{a = 1.691767707, d = 2.248911615}

(3)

plot(eval(sol, ad), x=0..eval(d, ad))

 

 

 

Download catenary_2.mw


foot note : with h1=50, h2=70 and L=160 fsolve can't find a solution (I killed it after one minute).
I suggest you to make the problem dimensionless to make it easier to solve.
For instance, let C > 0, and denote

(a, d) = solution of the problem y(0)=h0, y(d)=h1, length=L

then the solution of the problem

y(0)=C*h0, y(d)=C*h1, length=L*C

is 

(a/C, d*C)



Last remark: In a practical problem, a would be an input of the problem (a depends on the properties of the wire: lineic mass and Young modulus) and the length L of the cable would be an output, not an input.

3 4 5 6 7 8 9 Last Page 5 of 28