Carl Love

Carl Love

23928 Reputation

25 Badges

9 years, 309 days
Natick, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity

These are answers submitted by Carl Love

Like this:

z:= [
    0, 0, 0, 0, 0, 0, .689376362, 1.378752724, 2.068129087, 2.757505449, 0, 1.02920355,
    2.0584071, 3.087610649, 4.116814199, 0, 1.216469264, 2.432938529, 3.649407793,
    4.865877057, 0, 1.325720912, 2.651441823, 3.977162735, 5.302883646
    [ListTools:-LengthSplit](z,5), axis[1,2]= [tickmarks= ([$1..5]=~[$0..4])]


The Maple and Mathematica solutions that you show are mathematically equivalent. Both will give nonreal values for any real values of the a's and b's (unless a[1,1]*b[0,2] = a[0,2]*b[1,1], in which case the discriminant is 0). The expression under the square root sign cannot be positive for any real values of the parameters because it factors as (-1)*(a[1,1]*b[0,2] - a[0,2]*b[1,1])^2.

The absence of from an expression is not a reliable indicator that the expression is real-valued. The presence of is not a reliable indicator that it is not real-valued.

For rational-coefficient polynomial equations (such as yours) with parameters, you can get solve to break down all the relevant real intervals of the parameters which produce real solutions by adding the options parametric and real:

    k^2*a[1,1]^2 + k^2*b[1,1]^2 + 4*k*a[0,2]*a[1,1] + 
    4*k*b[0,2]*b[1,1] + 4*a[0,2]^2 + 4*b[0, 2]^2 = 0, 
    k, parametric, real

Note that the solution below covers all degenerate cases, i.e., when the equation is not really quadratic because the coefficient of k^2 is 0 and even when the coefficient of k is also 0. The parametric option always considers those cases.


I agree that a list or sequence of vectors would be the preferred format for programmatic usage. For a displayable table, I prefer using a matrix. (I especially like the way that matricies display with their entries centered in both dimensions.) The code below does both: My procedure DiffTab returns the "raw" list of vectors, the same as Tom Leslie's doDiff. The procedures PadV and DisplayTab turn those vectors into a neatly displayable matrix.

PadV:= proc(V::Vector, n::posint)
local k:= numelems(V), OP:= j-> n-k-1+2*j;
    Vector(2*n-1, applyop~(OP, 1, op(2,V)), fill= ``)
end proc
DisplayTab:= (DT::list(Vector), X::Vector)-> `<|>`(PadV~([X, DT[]], numelems(X))[])
DiffTab:= proc(Y::Vector)
local k, r:= Y, R:= table([1=r]);
    for k from 2 to numelems(Y) do r:= r[2..] - r[..-2];  R[k]:= r od;
    convert(R, list)
end proc
X:= <[$1..5]>:
Y:= 1/~X:
DT:= DiffTab(Y);
<<` X` | ` Y` | [``$4][]>, DisplayTab(DT,X)>;


You need to change xy to x*y.

Note the diffferent output from these two:

solve(y(1)^2 = 2., y(1));

1.414213562, -1.414213562

{ fsolve(y(1)^2 = 2., y(1)) };
             { }

The fsolve returns no output (which I emphasized by putting the command in { }), not even an error message. Technically, the fsolve output is called NULL.

An expression like y(1) is called in Maple an unevaluated function call (or sometimes just a "function"). So, you can use a function as a variable for solve but not for fsolve. So, your call T__v_x(1000) is producing NULL output in a position where some other command is expecting a number.


Do you mean something like this?

#x and y coordinates of grid:
GX:= [seq](0..3, 1/3):  GY:= [seq](1..3, 1/3)
#x and y coordinates of highlighted points:
PX:= [seq](1/4..11/4, 1/4):  PY:= [seq](5/4..11/4, 1/4)
        [[GX[1],GY[1]], [GX[1],GY[-1]], [GX[-1],GY[-1]], [GX[-1],GY[1]], [GX[1],GY[1]]],    
        color= COLOR(RGB, .3$3), thickness= 2
            seq([[GX[1],y],[GX[-1],y]], y= GY[2..-2]), 
            seq([[x,GY[1]],[x,GY[-1]]], x= GX[2..-2])
        color= COLOR(RGB, .5$3), thickness= 0
        [seq](seq([x,y], x= PX), y= PY),
        style= point, symbol= solidcircle, symbolsize= 8,
        color= red
    axes= none, scaling= constrained

Your general solution isn't general enough. If f is *any* differentiable function such that f(0)=1, then x0 = y0 = f is a solution to your equation. 

A basic fact about systems of equations is that the number of things that can be solved for is at most the number of equations. 

If is the list of expressions, then all that you need to do is

R:= Grid:-Map(int, L, x= 0..1);

For numeric integration, do 

R:= Grid:-Map(int, L, x= 0..1, numeric);

The command will automatically apply some rudimentary load balancing.

It's possible to get the exact solutions for all 6 functions by direct use (i.e., without using dsolve) of the Method of Undetermined Coefficients. If you try to do this with dsolve, it will get the first 5, but the 6th took longer than I cared to wait. I killed it after more than an hour. The method below gets it in seconds.

Direct application of Method of Undetermined Coefficients


#right sides of the 6 odes:
RHS:= [
    0, -v__0^3, -3*v__0^2*v__1, -3*v__0*v__1^2 - 3*v__0^2*v__2,
    -6*v__0*v__1*v__2 - 3*v__0^2*v__3 - v__1^3,
    -6*v__0*v__1*v__3 - 3*v__1^2*v__2 - 3*v__0^2*v__4 - 3*v__0*v__2^2







ICs:= [[1,0], [0,0]$5] #initial conditions of the 6 odes

for k to nops(RHS) do
    S:= C1*sin(t) + C2*cos(t); #general homogenous solution
    #Perform method of undetermined coefficients termwise:
    rs:= combine(expand(RHS[k]));
    for f in indets(rs, specfunc({sin,cos})) do #for each term on right side
        C:= coeff(rs, f);
        d:= degree(C,t);
        o:= op(f);
        P:= add(a[i]*t^i, i= 0..d)*sin(o) + add(b[i]*t^i, i= 0..d)*cos(o);
        if o=t then P*= t fi; #extra power of t   
        S+= eval(P, solve(identity(diff(P,t$2)+P=C*f, t)))
    #Apply initial conditions to determine C1 and C2:
    v__||(k-1):= (combine@expand@eval)(S, solve(eval([S, diff(S,t)], t=0)=~ICs[k]));
    print(v[k-1] = v__||(k-1))

v[0] = cos(t)

v[1] = -(1/32)*cos(t)-(3/8)*sin(t)*t+(1/32)*cos(3*t)

v[2] = (23/1024)*cos(t)+(3/32)*sin(t)*t-(3/128)*cos(3*t)+(1/1024)*cos(5*t)-(9/128)*cos(t)*t^2-(9/256)*t*sin(3*t)

v[3] = -(547/32768)*cos(t)+(9/1024)*sin(t)*t^3-(207/4096)*sin(t)*t+(135/4096)*cos(t)*t^2+(279/8192)*t*sin(3*t)-(81/4096)*t^2*cos(3*t)+(297/16384)*cos(3*t)-(3/2048)*cos(5*t)+(1/32768)*cos(7*t)-(15/8192)*t*sin(5*t)

v[4] = (1/1048576)*cos(9*t)-(99/16384)*sin(t)*t^3-(9/131072)*cos(7*t)+(1539/65536)*t^2*cos(3*t)+(825/262144)*t*sin(5*t)-(1359/65536)*cos(t)*t^2-(3915/131072)*t*sin(3*t)+(883/524288)*cos(5*t)-(15121/1048576)*cos(3*t)+(8997/262144)*sin(t)*t+(6713/524288)*cos(t)+(27/32768)*cos(t)*t^4-(225/131072)*t^2*cos(5*t)+(243/32768)*t^3*sin(3*t)-(21/262144)*t*sin(7*t)

v[5] = -(3/1048576)*cos(9*t)+(4635/1048576)*sin(t)*t^3-(81/1310720)*sin(t)*t^5+(1757/16777216)*cos(7*t)-(96795/4194304)*t^2*cos(3*t)-(16575/4194304)*t*sin(5*t)+(15777/1048576)*cos(t)*t^2+(108243/4194304)*t*sin(3*t)-(921/524288)*cos(5*t)+(394701/33554432)*cos(3*t)-(108357/4194304)*sin(t)*t-(42397/4194304)*cos(t)-(441/4194304)*t^2*cos(7*t)-(27/8388608)*t*sin(9*t)+(2187/1048576)*t^4*cos(3*t)+(1125/1048576)*t^3*sin(5*t)+(1/33554432)*cos(11*t)-(783/1048576)*cos(t)*t^4+(6975/2097152)*t^2*cos(5*t)-(10935/1048576)*t^3*sin(3*t)+(1659/8388608)*t*sin(7*t)



The plot command is plotting over a real interval. It can't be expected to find on its own every turn of a function, especially considering that there could be infinitely many of them (such as for the Weierstrass function[*1]). If there are values of n that you want to make sure are included, you can use the sample option. In this case, we might as well make the sample all the integers in the interval:

f:= n-> ceil(sqrt(4*floor(n))) - floor(sqrt(2*floor(n))) - 1:

Edit: My intention was to use the OP's original function, but I mistakenly copy and pasted Tom Leslie's modification of the original. So, make that

f:= n-> ceil(sqrt(4*n)) - floor(sqrt(2*n)) - 1:
plot(f, 10...100, sample= [$10..100]);

[*1] Weierstrass function: f:= x-> Sum(cos(3^n*x)/2^n, n= 0..infinity). It is continuous everywhere and differentiable nowhere.

You should use for derivatives in initial and boundary conditions. I don't know why that is, but it was the only change that I needed to make to your code. Sometimes the diff form works, but the form always works and is a lot less to type.

sys:= {
    diff(phi(eta), eta$2) + 5.261282735*f(eta)*diff(phi(eta), eta) - 
    2.630641368*phi(eta) = 0,
    1.059704409*diff(theta(eta), eta$2) + 6.176017503*f(eta)*diff(theta(eta), eta) 
    + 21.03607964*diff(f(eta), eta$2) + 0.5*phi(eta) = 0,
    diff(f(eta), eta$4) - 1.052256547*diff(f(eta), eta)*diff(f(eta), eta$2) + 
    1.052256547*f(eta)*diff(f(eta), eta$3) + 5.165076420*diff(theta(eta), eta) + 
    5.261282735*diff(phi(eta), eta) = 0,

    D(phi)(0) = 1 + 0.5*(D@@2)(f)(0),
    D(phi)(1) = 0.5*(D@@2)(f)(1), 
    f(0) = -0.5, f(1) = 0.5, 
    phi(0) = 1, phi(1) = 0, 
    theta(0) = 1, theta(1) = 0
dsol:= dsolve(sys, numeric):
    dsol, `[]`~(eta, [f, phi, theta](eta)), eta= 0..1, 
    legend= [f, phi, theta]

This is an attempt to Answer your followup Question.

A single summation with two indices, say, Sum(f(i, j), i+j= 0..N), is equivalent to the nested summation

Sum(Sum(f(i, j-i), i= 0..j), j= 0..N)

So, if you define your function as 

P:= (x,y)-> Sum(Sum(alpha[i, j-i]*x^i*y^(j-i), i= 0..j), j= 0..N)

then you can differentiate it like any usual function of x and y.

You began your worksheet by assigning numeric values to "symbols". Then you used those symbols in some formulas, and then you obtained or attempted to obtain some purely numeric solutions, such as plots. This is almost always a bad order to do things---whether you're using a symbolic language like Maple, a numeric language like Matlab, or just doing it with pen and paper: It may seem convenient at first, but it becomes inconvenient when you want to change something. And, if you're using a symbolic language, that order may prevent you from using some of the symbolic features. For example, in the worksheet below, I use some of the basic commands for manipulating symbolic polynomials. These don't work as well when the polynomials contain numeric constants with decimal points (called floats). (They still work, but not as well.)

The better approach is to start with a fully symbolic presentation of the problem and then incorporate the numerics at the end.

The central formula (equation) of your worksheet (the one in the solve command) can be easily converted to a polynomial in x. And, if you keep it symbolic, most above-average high-school algebra 2 students could do that conversion with pen and paper. If decimal numbers are used before this step, the result becomes a difficult-to-read mess. There are many benefits to using a polynomial. For example, fsolve automatically returns multiple solutions for polynomials, but not for other equations.

The worksheet below runs many, many times faster than the original approach of solving for 100 values of t.

NULLLCC Resonant Converter - Frequency Finder


Idea, technical formulas, and numeric data by


Maple code by Carl Love <>




#Very important: There are no decimal points
#used in the input *anywhere* in this worksheet!
params:= [
    #These are things that you set at the top of your worksheet. Note that I'm not
    #assigning these values *yet*, just declaring them and loosely associating
    #them with some symbols.
    V__out= 50,
    omega__line= 100*Pi,
    #t= ... Don't specify any t values at this time!
    #But we keep in mind that t goes from 0 to 1/100.
    theta= omega__line*t,
    V__line= 85*sqrt(2)*sin(theta),
    #m= ... I'm holding off on giving m a value here;
    #       I'll leave that up to you.
    Q__s= 16/5*sin(theta)^2,

    #I added this one because it helps simplify the equation:
    VR= V__line/V__out #"Voltage Ratio"
tr:= 0..1/100 #range of t.

#This is your original equation, except that I substituted my VR.
eq:= 1/VR = 1/sqrt((-m*x^2 + m + 1)^2 + (Q__s*(x - 1/x))^~2);

1/VR = 1/((-m*x^2+m+1)^2+Q__s^2*(x-1/x)^2)^(1/2)

#Convert it to a polynomial in x. Although that could be done by Maple, I think in this
#case it's easier to just use high-school algebra. The steps are
#    1. Reciprocate both sides.
#    2. Square both sides.
#    3. Subtract left side from right side (obtaining an expression implicitly
#       equated to 0).
#    4. Multiply by x^2.
#    5. Expand.
#    6. Collect powers of x.

#Those 6 steps are all done by this single line of code:
eqP:= collect(expand(x^2*((rhs-lhs)@map)(`^` , eq, -2)), x);


#Note that all powers of x are even. Since we're only interested in positive solutions,
#we can substitute X for x^2 *now*, then we'll take the square root of the numeric
#solutions *after* they're obtained.
eqP2:= simplify(eqP, {x^2= X});


#This "inner" procedure calls fsolve, not solve, and returns all
#solutions X >= 0. Note that m must be specified via an "outer" procedure.
FS_inner:= proc(t, X::name:= :-X)
option remember;
local r;
    if not t::realcons then return 'procname'(args) fi;
    r:= [fsolve](_P, X= 0..infinity, 'avoid'= {X=0});
    if not r::list(numeric) then
        printf("No solution found for t = %a.\n", t);
end proc

#The outer "wrapper" that handles m:
FS:= subs(
    _inner= eval(FS_inner),
    _P= subs(m= _m, eval[recurse](eqP2, params)),
    proc(m) option remember; rcurry(subs(_m= m, _inner), _rest) end proc    

#Operator for the final solutions: The square roots of the min and max of the
#polynomial solutions:
x_min_max:= m-> sqrt@~[min,max]@~FS(m, _rest):

#Plot for m=1:
    x_min_max(1), tr,
    title= "LCC Gain Frequency",
    labels= [t,x], labelfont= [TIMES, BOLD, 14],
    color= [red, green], thickness= 2, view= [tr, 0..2],
    caption= typeset(cat(` `$9, m) = 1),
    captionfont= [HELVETICA, BOLD, 16]




If you'd like the horizontal axis to be scaled by V__line (as VV did), I can make a small adjustment to do that.

You could use the top-level longstanding galois command for the information that you want.

galois(x^5 + 20*x + 32);
"5T2", {"5:2", "D(5)"}, "+", 10, {"(1 2 3 4 5)", "(1 4)(2 3)"}

The help page ?galois describes what all that means.

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

5 6 7 8 9 10 11 Last Page 7 of 352