Carl Love

Carl Love

28035 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

As Kitonum points out, this is done by the command GraphTheory:-SequenceGraph. However, this will return a graph with the vertices generically labeled from 1 to n. I got the impression from your Question that you wanted to use your own vertex labels. So here is a procedure to do that. The procedure takes two arguments. The first is the list of degrees. The second argument is optional and is the corresponding list of vertices. If omitted, generic 1-to-n labels will be used. The procedure returns FAIL if the degree sequence can't possibly be made into a graph.

RealizeDegreeSequence:= proc(
     Degrees::list(nonnegint),
     Vertices::list({symbol,string,integer}):= [$1..nops(Degrees)]
)
local
     Edges:= table(), e:= 0,
     D:= Vector(Degrees),
     V:= Vector(Vertices),
     P, k
;
     if max(D) >= op(1,D) then  return FAIL  end if;
     do
          (D,P):= sort(D, 'output'= ['sorted', 'permutation']);
          V:= V[P];
          if D[-1] = 0 then  break  end if;
          for k from 2 to D[-1]+1 do
               if D[-k] = 0 then  return FAIL  end if;
               D[-k]:= D[-k] - 1;
               e:= e+1;
               Edges[e]:= {V[-1],V[-k]}
          end do;
          D[-1]:= 0
     end do;
     GraphTheory:-Graph(Vertices, convert(Edges, set))
end proc:

I'd rather do the example with Sum rather than sum because the evaluation of the sums gives results that obsfucate the example, but I wrote the code so that it works for either. Using subsindets, it's a short one-liner.

 

S:= 2*(Sum(A*cos(omega*t)^2+cos(omega*t)*B*sin(omega*t)-cos(omega*t)*ymeas(t), t));

2*(Sum(A*cos(omega*t)^2+cos(omega*t)*B*sin(omega*t)-cos(omega*t)*ymeas(t), t))

subsindets[flat](S, specfunc({sum,Sum}), S-> map(op(0,S), op(S)));

2*(Sum(A*cos(omega*t)^2, t))+2*(Sum(cos(omega*t)*B*sin(omega*t), t))+2*(Sum(-cos(omega*t)*ymeas(t), t))

 

 

Download subsindets.mw

I found two substantial differences between your Mathematica code and your Maple code (from your PDF versions). In the Mathematica code, the initial condition for h is

G:= 0.005*exp(-(x-0.7)^2/(1/500)) + hin;

In the Maple code it's

F:= 0.01*exp(-(x-0.4)^2/(1/300)) + hin;

The other difference is that in the Mathematica code, the upper bound of is w = 1, whereas in the Maple it's l = 2.

If I use your Mathematica conditions in the Maple code, the animation runs much faster.

 

I redid your worksheet in a much cleaner style (even though it was already quite well organized), which is much less prone to user errors. All floating-point constants are in a list of parameters that is not applied to the system until needed. Thus you can see the PDEs in their fully symbolic form. I changed zita to zeta[0] and gama to gamma (and made gamma local) so that these would print nicer. 

Then I used numeric differentiation (fdiff) applied to the pdsolve solution to compute the residuals of the two PDEs over a grid of 18 interior points. The residuals seem reasonably small to me, except near (x,t) = (1,5). Also, read the help page ?pdsolve,numeric,errorcontrol.

 

restart:

 

Digits:=15:
interface(rtablesize= 30):

local gamma:

 

PDEtools:-declare(h(x,t),u(x,t)):

params:= [
     g= 9.81,
     beta= 0.65,
     L= 0.001,
     zeta[0]= evalf(35.1*Pi/180),
     zeta[1]= evalf(32.9*Pi/180),
     zeta[2]= evalf(42*Pi/180),
     mu1= tan(zeta[1]),
     mu2= tan(zeta[2]),
     gamma= (mu2-mu1)/(tan(zeta[0])-mu1) - 1,
     B= beta*sqrt(g*cos(zeta[0]))/(L*gamma),
     nu= 2*gamma*L*sqrt(g)*sin(zeta[0])/(9*beta*sqrt(cos(zeta[0]))),
     hstop= L*gamma,
     hin= 1.2*hstop,

     uin= B*hin^(3/2),
     l= 2 #upperbound of x
]:

'params' = <params[]>;
 

mu:= mu1+(mu2-mu1)/(1+beta*h(x,t)^(3/2)*sqrt(g*cos(zeta[0]))/(L*u(x,t)));

 

 

h(x, t)*`will now be displayed as`*h

u(x, t)*`will now be displayed as`*u

params = (Matrix(15, 1, {(1, 1) = g = 9.81, (2, 1) = beta = .65, (3, 1) = L = 0.1e-2, (4, 1) = Zeta[0] = .612610567450009, (5, 1) = Zeta[1] = .574213323906135, (6, 1) = Zeta[2] = .733038285837617, (7, 1) = mu1 = tan(Zeta[1]), (8, 1) = mu2 = tan(Zeta[2]), (9, 1) = gamma = (mu2-mu1)/(tan(Zeta[0])-mu1)-1, (10, 1) = B = beta*sqrt(g*cos(Zeta[0]))/(gamma*L), (11, 1) = nu = (2/9)*gamma*L*sqrt(g)*sin(Zeta[0])/(beta*sqrt(cos(Zeta[0]))), (12, 1) = hstop = gamma*L, (13, 1) = hin = 1.2*hstop, (14, 1) = uin = B*hin^(3/2), (15, 1) = l = 2}))

mu1+(mu2-mu1)/(1+beta*h(x, t)^(3/2)*(g*cos(zeta[0]))^(1/2)/(L*u(x, t)))

F:= 0.01*exp(-(x-0.4)^2/(1/300))+hin;

plot(eval[recurse](F, params), x= 0..2);

 

0.1e-1*exp(-300*(x-.4)^2)+hin

pde1:= Diff(h(x,t),t)+Diff(h(x,t)*u(x,t),x)=0;

Diff(h(x, t), t)+Diff(h(x, t)*u(x, t), x) = 0

pde2:=
     Diff(h(x,t)*u(x,t), t) + Diff(h(x,t)*u(x,t)^2, x) =
     h(x,t)*g*sin(zeta[0]) -
     mu*h(x,t)*g*cos(zeta[0]) -
     Diff((1/2)*(h(x,t))^2*g*cos(zeta[0]),x) +
     Diff(nu*h(x,t)^(3/2)*Diff(u(x,t),x),x)
;

 

Diff(h(x, t)*u(x, t), t)+Diff(h(x, t)*u(x, t)^2, x) = h(x, t)*g*sin(zeta[0])-(mu1+(mu2-mu1)/(1+beta*h(x, t)^(3/2)*(g*cos(zeta[0]))^(1/2)/(L*u(x, t))))*h(x, t)*g*cos(zeta[0])-(Diff((1/2)*h(x, t)^2*g*cos(zeta[0]), x))+Diff(nu*h(x, t)^(3/2)*(Diff(u(x, t), x)), x)

sys:= [pde1,pde2]:

ibc:= [
     h(x,0)=F, h(0,t)=h(l,t), u(x,0)=uin,
     D[1](u)(0,t) = D[1](u)(l,t), u(0,t)=u(l,t)
]:

'ibc' = <ibc[]>;

 

ibc = (Matrix(5, 1, {(1, 1) = h(x, 0) = 0.1e-1*exp(-300*(x-.4)^2)+hin, (2, 1) = h(0, t) = h(l, t), (3, 1) = u(x, 0) = uin, (4, 1) = (D[1](u))(0, t) = (D[1](u))(l, t), (5, 1) = u(0, t) = u(l, t)}))

sol:= pdsolve(
     eval[recurse](value([sys,ibc]), params)[], [h,u], numeric,
     time= t, spacestep= 0.001, timestep= 0.06,
     range= eval(0..l, params),
     compile= true,
     optimize= true
);

 

module () local INFO; export plot, plot3d, animate, value, settings; option `Copyright (c) 2001 by Waterloo Maple Inc. All rights reserved.`; end module

sol:-animate(
     h(x,t), t= 0..6.2, thickness= 2, frames= 100,
     color= blue, view= [default,0..0.1]
);

 

(U,H):= eval([u(x,t),h(x,t)], sol:-value(output= listprocedure))[];

`[Length of output exceeds limit of 1000000]`

Residuals:= unapply(
     subs(
          diff(u(x,t),x$2)= fdiff(U, [1,1], [x,t]),
          diff(u(x,t),x)= fdiff(U, [1], [x,t]),
          diff(u(x,t),t)= fdiff(U, [2], [x,t]),
          diff(h(x,t),x)= fdiff(H, [1], [x,t]),
          diff(h(x,t),t)= fdiff(H, [2], [x,t]),
          u(x,t)= U(x,t), h(x,t)= H(x,t),
          (lhs-rhs)~(eval[recurse](value(sys), params))
     ), [x,t],
     proc_options= remember
);     

proc (x, t) option remember; [fdiff(H, [2], [x, t])+fdiff(H, [1], [x, t])*U(x, t)+H(x, t)*fdiff(U, [1], [x, t]), fdiff(H, [2], [x, t])*U(x, t)+H(x, t)*fdiff(U, [2], [x, t])+fdiff(H, [1], [x, t])*U(x, t)^2+2*H(x, t)*U(x, t)*fdiff(U, [1], [x, t])-5.64080152254456*H(x, t)+8.02604872793949*(.646929012633441+.253475031664396/(1+1841.46832379882*H(x, t)^(3/2)/U(x, t)))*H(x, t)+8.02604872793949*H(x, t)*fdiff(H, [1], [x, t])-0.361033442303963e-2*H(x, t)^(1/2)*fdiff(U, [1], [x, t])*fdiff(H, [1], [x, t])-0.240688961535974e-2*H(x, t)^(3/2)*fdiff(U, [1, 1], [x, t])] end proc

for X from .5 by .5 to 1.5 do
     for T from 1 by 1 to 6 do
           Residuals(X,T)
     end do
end do;

op(4, eval(Residuals));

table( [( 1.5, 1 ) = [HFloat(-8.7e-11), HFloat(-1.2522802239622877e-11)], ( .5, 5 ) = [HFloat(2.85953750500176e-5), HFloat(-0.001490879050560265)], ( 1.0, 1 ) = [HFloat(1.932828276029985e-6), HFloat(4.135527008128567e-5)], ( .5, 1 ) = [HFloat(2.287801319130747e-6), HFloat(-2.0692549446668662e-5)], ( 1.5, 5 ) = [HFloat(2.053365866274521e-6), HFloat(1.0652653516122965e-4)], ( .5, 3 ) = [HFloat(-1.1647776369834434e-7), HFloat(1.318382691491558e-7)], ( .5, 4 ) = [HFloat(6.435492967658586e-4), HFloat(-0.012572282194018345)], ( 1.0, 2 ) = [HFloat(-6.725756623120064e-4), HFloat(5.915982663221885e-4)], ( 1.0, 4 ) = [HFloat(7.606652077449086e-7), HFloat(-3.3512252416319854e-5)], ( 1.0, 3 ) = [HFloat(5.244656243204955e-5), HFloat(8.974099297601484e-5)], ( 1.5, 3 ) = [HFloat(0.016529087254550233), HFloat(0.02796219302039396)], ( .5, 6 ) = [HFloat(-8.402923824848941e-6), HFloat(1.1685994412418874e-4)], ( 1.5, 2 ) = [HFloat(1.2992914585305754e-5), HFloat(2.5073699512041575e-5)], ( .5, 2 ) = [HFloat(-5.696111396287452e-7), HFloat(1.7360380724026877e-6)], ( 1.0, 5 ) = [HFloat(0.19337670235490467), HFloat(0.06306079458023821)], ( 1.0, 6 ) = [HFloat(-4.413769663749155e-5), HFloat(-0.0015418020395394545)], ( 1.5, 4 ) = [HFloat(2.88274342677532e-5), HFloat(-1.0790727271419618e-4)], ( 1.5, 6 ) = [HFloat(4.155903691802692e-4), HFloat(-0.0021539603521344894)] ] )

 

 

Download PDE_residuals.mw

The Maple command for (infinite) series is sum:

sum(f, n= 2..infinity);

Maple returns infinity, meaning that the series diverges. To restrict to odd terms, do

sum(eval(f, n= 2*k+1), k= 1..infinity);

If you actually mean the convergence or divergence of the sequence rather than the series, use the limit command:

limit(f, n= infinity);

which returns 0.

Considering the triviality of this problem, I suggest that you look up the Limit Comparison Test in your textbook or on the Web. If you knew this test, you would be able to look at this problem and instantly know that the series diverges.

This Answer is to answer the OP's followup question requesting a procedure that determines whether an arbitrary polynomial can be produced as a linear combination of a set of polynomials, and, if so, returning the coefficients.

I decided that it would be better to just write a separate procedure for this rather than rolling it into my procedure LinearCombo from an earlier thread. So I wrote a little procedure PolyLinearCombo that calls my version of Kitonum's CoeffsOfLinComb. I changed the declaration-of-variables argument to be optional; by default, all nonconstant names will be considered variables. If you want to use parameters, then include the variables argument as a set.

For various reasons, it's better that a procedure always return the same number of return values. So, this procedure returns true and a list of coefficients if the second argument can be produced as a linear combination of the elements of the first argument, and it returns false, [] if there's no such combination.

I decided to disallow 0 as the second argument. If you want nontrivial coefficients to produce 0, then just call CoeffsOfLinComb.

CoeffsOfLinComb:= proc(
     L::list(polynom),
     V::set(name):= indets(L, And(name, Not(constant)))
)
local
     c, k, C:= {c[k] $ k= 1..nops(L)},
     S:= solve({coeffs(expand(`+`((C*~L)[])), V)}, C),
     F:= indets(rhs~(S), name) intersect C=~ 1
;
     eval([C[]], eval(S,F) union F)
end proc:


PolyLinearCombo:= proc(
     F::list(polynom),
     f::And(polynom, Not(identical(0))),
     V::set(name):= indets([F,f], And(name, Not(constant)))
)
local C:= CoeffsOfLinComb([f, F[]], V);
     if C[1]=0 then (false, []) else (true, -C[2..] /~ C[1]) end if
end proc:

Examples of use:

Dependent example:

F:= ['randpoly([x,y])' $ 22]:
f:= randpoly([x,y]):
C:= PolyLinearCombo(F, f, {x,y});

Independent example:

F:= ['randpoly([x,y])' $ 22]:
f:= randpoly([z,y]):
C:= PolyLinearCombo(F, f, {x,y,z});

The above example becomes dependent if z is treated as a parameter:

F:= ['randpoly([x,y])' $ 22]:
f:= randpoly([z,y]):
C:= PolyLinearCombo(F, f, {x,y});

I'd condsider this a bug: The value of Digits and other environment variables (which have been modified from their defaults) should be transcribed to all nodes. Another workaround is to put

Digits:= 30;

into your initialization file.

I agree with Roman that you should run the evalf also in parallel; timewise, it is a significant part of the computation. But I do not agree that you should necessarily use numeric integration. So, change the Mapped procedure from int to evalf@Int or evalf@int depending on whether you want to use numeric or symbolic integration.

Your problem is that pdf1(i), etc., will evaluate to a numeric value even if i is symbolic (i.e., even if i has no numeric value). There are two ways to correct this. The preferred way is to make the following the first executable line of pdf1:

if not x::numeric then return 'procname'(args) end if;

The crude way is to put unevaluation quotes around the call to pdf1:

Sum('pdf1'(i), i= 3..3),

etc.

There is no bug.

The subtleties of unevauation-quote syntax can be mind-boggling. In this case, I find it easiest to take advantage of the fact that Explore's first argument is automatically unevaluated (see showstat(Explore, 1) or op([1,1], eval(Explore:-ModuleApply))) by wrapping all plots into the call to Explore. This can be done in a one-liner. I split the following command over a few lines so that it would display nicely on MaplePrimes, but it'll easily fit on one line in a Maple worksheet:

Explore(
     plots:-display(`<|>`(plot~([a*x^2, a*x^2+1], x= -1..1, y= -3..3)[])),
     parameters= [a= -1.0..1.0]
);

Syntax highlights:

`<|>`(1,2) is equivalent to < 1 | 2 >; it creates a row vector, which is sufficient for a side-by-side array of plots.

L[] converts list into its corresponding sequence. So `<|>`([1,2][]) is equivalent to < 1 | 2 >.

f~([a,b], c) is equivalent to map(f, [a,b], c) is equivalent to [f(a,c), f(b,c)].

For set complements, you could overload the exponentiation operator so that a set with superscript would represent the set complement. Indeed, it could all be put into a module constructed so that all the operations are done with respect to a particular universal set. So, below is essentially the same Answer that I gave before, but all the notations are shorter. I made `+` represent set union and `*` for intersection. I no longer feel that this is too trivial---the overloaded operators a were the missing piece.


restart:

SetOperations:= proc(U::set)
     module()
     option package;
     export
          `^`:= proc(A::set, e::identical(c), $)
          option overload;
               U minus A
          end proc,
          
          Select:= (f::procedure)-> select(f, U),
          
          `+`:= proc(A::set, B::set)
          option overload;
               A union B
          end proc,

          `*`:= proc(A::set, B::set)
          option overload;
               A intersect B
          end proc,
       
          P:= proc(A::set, B::set:= U, $)
               nops(A intersect B)/nops(B)
          end proc          
     ;
     end module
end proc:

TwoDice:= SetOperations({combinat:-permute([$1..6, $1..6], 2)[]}):
with(TwoDice);

[`*`, `+`, P, Select, `^`]

H:= Select(x-> :-`+`(x[]) = 5);

{[1, 4], [2, 3], [3, 2], [4, 1]}

G:= Select(x-> x[2] < 3);

{[1, 1], [1, 2], [2, 1], [2, 2], [3, 1], [3, 2], [4, 1], [4, 2], [5, 1], [5, 2], [6, 1], [6, 2]}

P(H);

1/9

P(G);

1/3

P(G^c);

2/3

P(H+G);

7/18

P(H*G);

1/18

P(G,H);

1/2

P(H,G);

1/6

#Law of total probability illustration

P(H) = P(H,G)*P(G) + P(H,G^c)*P(G^c);

1/9 = 1/9


Download prob.mw

 

 

You could use eval statements rather than assignment statements, like this

restart:
Z:= Matrix(2,2,symbol= z):
for k to 2 do
     Z||k:= eval(Z, z[k,k]= 1)
end do:
Z, Z1, Z2;

Thus, the values of and z never change. You're only making "working copies" of Z

Use P[A]:= 0.3, etc. For conditionals, use P[A,B].

The syntax f1(p,b):= ... is not the correct way to define a function in Maple. The correct syntax is 

f1:= (p,b)-> 1/(p^2 + b^2)^2;

Then enter the integral as 

T1pbm:= int(q*f1(q,b)*ln((p+q)^2 +m^2), q= 0..infinity) assuming real;

Notice that int is with a lowercase i. The integral will take 5-10 seconds. Then do the same for the second integral. It will take several minutes.

Since the integrals will be done, the calls to Parts will be superfluous. Since your variables have no numeric values, the calls to evalf are superfluous. Finally, the constant Pi is spelled with a capital P in Maple; pi is just another variable.

The equals sign is not the assignment operator in Maple, so your line in the for loop is not setting the values of A[0] and A[1]. The assignment operator is colon-equals. So, you simply need to change the line in the for loop to 

A[n] := int(...)

Here's a completely different algorithm. By evaluating the equation

f - add(c[k]*g||k, k= 1..m) = 0  (*)

at different exact random numeric evaluations of the variables xy, ..., we obtain a set of constant-coefficient linear equations in unknowns c[1], ..., c[m]. This system is solved. If there's no solution, then there's no linear combination (guaranteed). If there's a solution, it's used in a symbolic attempt to reduce (*) to 0=0. If this is possible, then it's a linear combination (guaranteed).

This algorithm works over a class of expressions much larger than polynomials, and it handles transcendental coefficients. Here it is:

LinearCombo:= proc(B::list(algebraic), T::algebraic)
local
     V:= indets([B,T], And(name, Not(constant))),
     nV:= nops(V),
     c, #linear coeffiecients
     k, K, #indices
     nB:= nops(B),
     Eq:= add(c[k]*B[k], k= 1..nB) - T,
     Eqs,
     S, Z, R
;
     for K to 2 do
          Eqs:= {'eval(Eq, V=~ {'rand()' $ nV})' $ nB};
          if nops(Eqs) = nB then break end if;
     end do;
     if K = 3 then
          error "Almost certainly, all expressions are constant"
     end if;
     S:= solve(Eqs);
     if S=[][] then return false end if;
     Z:= indets(rhs~(S), And(name, Not(constant)))=~ 0;
     S:= eval(S,Z) union Z;
     R:= is(eval(Eq,S) = 0);
    `if`(R, eval([c[k] $ k= 1..nB], S), R)    
end proc:

The procedure is called liked this:

LinearCombo([y^2-y, y+x^2], x^2+y^2);

Other than errors for syntactically incorrect input or system problems, it returns one of four things:

  1. a list of coefficients that make a linear combination
  2. false, indicating that no linear combination exists
  3. rarely, FAIL, indicating the inability of is to simplify the coefficients
  4. the error message Almost certainly, all expressions are constant.

A positive result (1) is guaranteed correct (it's been verified with is). A false result can be wrong with a small probability if there is a highly convoluted relationship among coefficients containing transcendentals that is incorrectly labelled false by is. The result FAIL, of course, can never be wrong, and it represents a failure of is (see ?is). If the error "Almost certainly, all expressions are constant" is returned, there is an infintesimal probability that two sets of m random evaluations of the variables just happened to pick values such that an equation was redundant.

Comments welcomed, especially attempts to break the program or interesting experiences with transcendental functions, transcendental coefficients, trigonometric and other identities.

First 257 258 259 260 261 262 263 Last Page 259 of 395