Carl Love

Carl Love

28055 Reputation

25 Badges

12 years, 353 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

On the other hand, any place where a space is allowed, you can use multiple spaces, or line breaks, or any combination of those.

The command convert(..., binary) is not very useful for doing further arithmetic; it's mainly to display as output a binary representation of a number. An extremely efficient alternative is Bits:-Split, which returns a list of 0s and 1s. It's trivial to add a list with add (after Maple 2018, I believe) or `+`@op (in any Maple). To create your histogram, you could do

Statistics:-Histogram((`+`@op@Bits:-Split)~([$0..2^16]), binwidth= 1);

Another option that returns the same list is convert(..., base, 2). This is much slower than Bits:-Split (which is compiled code), but it's a command worth knowing because the 2 can be replaced by any positive integer other than 1; whereas Bits:-Split can only work with bases that are powers of 2. And it's not a very slow command; it's just relatively slow compared to the super-efficient Bits:-Split. The above command takes 2 seconds on my computer. If I use convert(..., base, 2), then it takes 7 seconds.

Both of these commands return the digits in order from least significant on the left to most significant on the right. Although that looks backward to us due to the way we usually write numbers, it is the form most convenient for further processing. It is only due to cultural rather than mathematical reasons that it looks weird. 

The dsolve has returned an equation whose left side is y(x) and whose right side is some expression containing but not y. We call that an explicit solution. You don't want to carry around that left side for the remaining computations. Change the line

C := eval(N__1, x = 200)

to

C:= eval[recurse](y(x), [N__1, x= 200])

There's no need to "estimate" the distance; it can be done exactly. Modify your distance function, like this:

dist:= proc(p::[integer,integer], q::[integer,integer])
local H;
    min((add@abs~)~([seq]([p[2]-H, p[1]-q[1], H-q[2]], H= [1,10])))
end proc  
:
Edgs:= {seq}([({op}@lhs~)(e), dist(rhs~(e)[])], e= combinat:-choose(pts2, 2)):
G1:= GraphTheory:-Graph(Edgs):
(distance, Tour):= GraphTheory:-TravelingSalesman(G1);
            distance, Tour := 38, [A, D, C, B, E, A]

I'm working on a method to highlight the trail with arrows.

Your troubles were caused by a hard-to-see typo. The subexpression beta(sin(alpha)/sin(beta)) occurs in the denominator of eq1. That beta is not multiplied by its following expression. To get implied multiplication in 2D Input, you need a space after beta. Or, you can use an explicit multiplication symbol *. If I make this correction and specify variable ranges to fsolve, then I get a solution close to (but not identical to) the one you mentioned earlier.

restart:

Digits:= 15:

params:= [
    Rrt= R/t, Ea= E/(1-nu^2),
    R=1511.6, t= 25.4, E=2e11, nu=0.29, F= 221e6
]:

#Abbreviations for common subexpressions:
macro(
    a= alpha, b= beta,
    S = sin(alpha)/sin(beta),
    PB= (3*Pi/2/beta)^2 - 1,
    T= tan(alpha - beta)
):
# Using the abbreviations reduces the chance of typos when entering the equations.
#

eq:= [
    sqrt(PB*(Pi - a + b*S^2)/12/S/(a - Pi/2500 - b*S/(1 + T^2/4)))/S - Rrt,
    PB/S^3 - 12*P*Rrt^3/Ea,
    1000*((1 - 1/S)/2/Rrt + P*Rrt*S/Ea + 4*P*Rrt^2*b*S^2*T/Pi/Ea - F/Ea)
]:
print~(eq): #for neater display
#
# Expressions are implicitly equated to 0.
#

(1/6)*3^(1/2)*(((9/4)*Pi^2/beta^2-1)*(beta*sin(alpha)^2/sin(beta)^2+Pi-alpha)*sin(beta)/(sin(alpha)*(alpha-(1/2500)*Pi-beta*sin(alpha)/(sin(beta)*(1+(1/4)*tan(alpha-beta)^2)))))^(1/2)*sin(beta)/sin(alpha)-Rrt

((9/4)*Pi^2/beta^2-1)*sin(beta)^3/sin(alpha)^3-12*P*Rrt^3/Ea

500*(1-sin(beta)/sin(alpha))/Rrt+1000*P*Rrt*sin(alpha)/(sin(beta)*Ea)+4000*P*Rrt^2*beta*sin(alpha)^2*tan(alpha-beta)/(sin(beta)^2*Pi*Ea)-1000*F/Ea

Eq:= eval[recurse](eq, params):

sol:= fsolve(Eq, {a= 0..Pi/2, b= 0..Pi/2, P= 0..infinity});

{P = 2065112.19495371, alpha = .944098810392359, beta = .934530721839835}

evalf(eval(Eq, sol)); #Check residuals.

[-0.217e-10, -0.2e-12, 0.58e-13]

[a,b]=~ evalf(eval([a,b]*180/Pi, sol)); #Convert angles to degrees.

[alpha = 54.0928772788039, beta = 53.5446661867369]

# Check Jacobian:
J:= Matrix((3,3), (i,j)-> evalf(eval(diff(Eq[i], [a,b,P][j]), sol)));

Matrix(3, 3, {(1, 1) = -5492.36319652914, (1, 2) = 5312.58704415908, (1, 3) = 0., (2, 1) = -51.9581308624461, (2, 2) = -.2740712344686, (2, 3) = -0.115826733650127e-4, (3, 1) = 47.4307787914329, (3, 2) = -47.1591891411326, (3, 3) = 0.4617082147e-6})

LinearAlgebra:-SingularValues(J);

Vector(3, {(1) = 7641.694050783875, (2) = 36.331927336981565, (3) = 0.7550434006e-6})

 

Download CriticalPressureFsolve.mw

If you

  1. use 15 Digits of precision,
  2. use fsolve instead of solve,
  3. use Pi instead of 3.1415, and
  4. search for complex solutions,

then the imaginary parts of those solutions will be very small, and you can ignore them. After you've entered the three equations, do this:

Digits:= 15:
PI:= Pi:
solC:= fsolve((Eq:= {eq||(1..3)}), complex);
         /                        8                      -11    
solC := { P = -1.44914053377011 10  + 1.81361988528669 10    I, 
         \                                                      
                                                 -21    
  alpha = -3.21040506122157 - 2.69567382732802 10    I, 

                                                -21  \ 
  beta = -2.50058386568901 - 1.52201908062012 10    I }
                                                     / 
solR:= eval(solC, I= 0); #Remove imaginary parts
         /                        8                             
solR := { P = -1.44914053377011 10 , alpha = -3.21040506122157, 
         \                                                      

                          \ 
  beta = -2.50058386568901 }
                          / 
evalf(eval(Eq, solR)); #Check residuals
       /       -12              -12             -10     \ 
      { -6.9 10    = 0., -3.3 10    = 0., 1.1 10    = 0. }
       \                                                / 

Red text is your input; blue text is Maple's output..

The command mtaylor does it. The code below assumes that y__i is some analytic expression in the variables x[1], x[2], x[3].

V:= [x[1],x[2],x[3]]:
ord:= 4:
sort(mtaylor(y__i, V, ord+1), V, ascending);

 

Like this:

V2:= -2*x + 4.95:
l__1:= 1.665:
plot(eval(V2, x= x - l__1), x= 
whatever);

Like this:

A:= {{a, b, c}, {a, b, d}};
add(add(x), x= A)

This works regardless of the size of A or the sizes of its elements. 

The easiest way to get all 22 real roots is

Student:-Calculus1:-Roots(cos(x) = x^2/1000);

If you use the decimal float coefficient 0.001, then you need to specify a finite interval:

Student:-Calculus1:-Roots(cos(x) = 0.001*x^2, x= -35..35);

I agree with Preben that DirectSearch:-SolveEquations is functioning as expected and as documented on its help page. You need to examine the residuals column, or set an appropriate option.

I suspect that the intended main function is exp(x)*cos(x) rather than exp(x*cos(x)). The integral for the intended area is then

int(abs(exp(x)*cos(x)), x= -Pi//2..Pi/2);

The abs (absolute value) makes no difference in this particular case, but it should be included in "area bounded by" problems to account for integrands that could be negative over some part of the interval of integration.

Does this do what you want?

L1:= [[seq1[i] $ i= 1..7]$24];
L2:= [$1..24]:
L3:= `~`[`+`]~(L1, [0, (L2[2..]-~1)[]]);

There is a simple bug in `evalf/doublefactorial`(z) only for z=1. The bug is in line 5 of showstat(`evalf/doublefactorial`): The 1.+z should be 1.. (And I wonder why whoever wrote this thought that z=1 needed to be treated as a special case. It seems to me that only z=0 needs special treatment.)

You can work around the bug by setting 

`evalf/doublefactorial`(1.):= 1.:

at the start. This is called a remember-table assignment.

This bug only affects numerical evaluation (via using evalf or floats). Symbolics aren't affected.

These three things that you believe are either not true or not valuable, and if you did accomplish them (which I could easily tell you how to do), you'd find out after that they aren't as good as the better ways that I show below:

  1. "But since the list contains the values as an expression, for example DQV[0, 0, 0, 0] = 0. In such cases, I cannot perform mathematical operations for the solution points results."
  2. "if I can delete the DQV[0,0,0,0] from the expression mentioned above, I'll end with the numeric values list, which can then be converted into an array and make things easier for me."
  3. "I'm struggling to delete the DQV part of the expression in one go (manually deletion for each expression, it's very time consuming to delete for 1500 variables)."

The first isn't true because you can easily use eval(expr, S) (where expr is any expression containing any subset of the variables and S is defined below) to numerically evaluate any mathematical expression of the variables. An example of using eval is given in example 3 below.

The second isn't valuable because it can easily be done without deleting anything, as shown in examples 2 and 4 below.

The third isn't valuable. In other words, there's no need to manually edit any output. I can't imagine anything useful that could be done with an LPSolve solution that couldn't be done with some simple commands (usually just 1 line) and no manual editing.

#Extract the list of variable = value equations:
S:= A:-Results("solutionpoint"):

#Example 1: Find all variables that are equal to 0:
ZeroVs:= lhs~(select(e-> rhs(e)=0, S));

#Example 2: Create a table (essentially an array) of all DQV variables
#(You shouldn't name it DQV, but Dqv is fine):
Dqv:= table(((op@lhs)=rhs)~(select(e-> type(lhs(e), specindex(DQV)), S))):

Dqv[1,1,1,0]; #etc.

#Example 3: Let C be your set of constraints. The following finds the subset that are
#inequality constraints that are at their limit at the current solution:
ActiveCs:= select(c-> c::`<=` and eval((lhs-rhs)(c), S) = 0, C);

#Example 4: Like example 2, but create a table of tables of all variables:
All:= table():
for e in S do
    if lhs(e)::indexed then All[op(0,lhs(e)][op(lhs(e))]:= rhs(e)
    else All[lhs(e)]:= rhs(e)
    fi
od:

All[DQV][1,1,1,0]; #etc.

I'd be happy to give more such examples if you say more precisely (mathematically) what you want. And you should ask here if you ever feel a need to manually edit any Maple output to produce input for other Maple commands.

Change the for line to

for ind, val in eval(F_vRk) do

Tables, procedures, and some modules have a property called last name evaluation. See ?last_name_eval.

First 49 50 51 52 53 54 55 Last Page 51 of 395