## dharr

Dr. David Harrington

## 6292 Reputation

19 years, 347 days
University of Victoria
Professor or university staff

## Social Networks and Content at Maplesoft.com

I am a professor of chemistry at the University of Victoria, BC, Canada, where my research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

## File restriction...

@Ronan My survey of the various filetypes used (.off, .ply, .stl etc) is that they are all for 3D objects and do not contain other information, with the possible exception of polygon face "textures" and "materials", which are for gaming applications and 3D printing respectively. So I think this is not a Maple issue. The face coloring I am less certain about, and the ability to support color is different for different files. But Maple doesn't seem to attempt to export color for any of them, which may be just that most 3D objects have a single color and not separate colors for individual faces. That can be fixed by messing with the PLOT3D structure, but that is significant effort.

## modified @sursumCorda...

@dharr @acer's routine is brilliant for a single sum, but does not take advantage of efficiencies if you need a sequence of sums. @sursumCorda's approach can be adapted to give an efficient way to generate a sequence.

 > restart;

Modification of @sursumCorda's routine that works out a term for each permutation of a partition.

This version  (i) accumulates the sums and partitions as we go (`combstruct/allstructs/Partition/allpartk` uses option remember), and

(ii) does not work out all the permutations, just takes one term for each partition and multiplies by the number of permutations

 > sums := proc(m::posint)   local n, s, parts, part, ss;   ss:=Vector(m);           s := 0:   for n to m do     parts := `combstruct/allstructs/Partition/allpartk`(n,3); # partitions of n with 3 parts     s := s + add(combinat:-numbperm(part)/mul(part), part in parts);     ss[n] := s;   end do;   ss end proc:
 > L:=CodeTools:-Usage(sums(200),iterations=4):

memory used=357.30MiB, alloc change=72.00MiB, cpu time=406.25ms, real time=3.39s, gc time=23.44ms

Further efficiencies are possible by generating information for each partition iteratively only once..
The partitions themselves do not have to be generated.
If  are the partitions of  with  parts, then these may be found by

(i) appending 1 to each of the

The product of the parts is unchanged.
The number of 1's increases by 1
The number of permutations is multiplied by , where  is the number of 1's in the old partition,
and  is the number of parts in the old partition.
For  quantity can be calculated as  for

and

(ii) adding 1 to each part in each of the

The new product is , with  for  and  for
The new  is
The number of 1's is now zero.

The number of permutations is unchanged.

 > sums2:=proc(m::posint)         description "outputs length m Vector with V[n] = sum(1/(i*j*k), 1 <= (i,j,k) <= n)";         local prods, m1, nperms, p2, n, s, ss, key, nn, nm1, nm2, nm3;         if m = 1 then return Vector([0]) end if;         if m = 2 then return Vector([0, 0]) end if;         if m = 3 then return Vector([0, 0, 1]) end if;         # initialize Array         # prods = products of partition parts, e.g. (2,2) for 2 in 2 parts = [1,1]; prod =1         prods := Array(1..4, 2..3, {                 (1,2) = Vector(),    (1,3) = Vector(),                 (2,2) = Vector([1]), (2,3) = Vector(),                 (3,2) = Vector([2]), (3,3) = Vector([1])});         # multiplicity of 1's, e.g., (2,2) has 2 1's         m1 := Array(1..4, 2..3, {                 (1,2) = Vector(),    (1,3) = Vector(),                 (2,2) = Vector([2]), (2,3) = Vector(),                 (3,2) = Vector([1]), (3,3) = Vector([3])});         # number of permutations, e.g. 2 ways for (3,2) = [1,2] or [2,1]         nperms := Array(1..4, 2..3, {                 (1,2) = Vector(),    (1,3) = Vector(),                 (2,2) = Vector([1]), (2,3) = Vector(),                 (3,2) = Vector([2]), (3,3) = Vector([1])});         # for partitions with 3 parts only, p2 is sum of (products of two parts)         # for [a,b,c], p2 = a*b + b*c +a*c                 p2 := Array(1..4, {(1) =  Vector(), (2) =  Vector(), (3) =Vector([3])});         s := 1;         ss := Vector(m, {1 = 0, 2 = 0, 3 = s});         # key has pointers to Array rows; reindexing allows reuse of rows in rotationg fashion         key := [1, 2, 3, 4];                 for n from 4 to m do                 nm3, nm2, nm1, nn := key[]; # mn3 = n-3 etc                 # update Array for this row                 # first part of Vector is for appending 1 to p(n-1,k-1)                 # second part of Vector is for adding 1's to p(n-k,k)                 # k = 2;                         # (n-1,1) = [n-1] -> [n-1,1]; adding 1's -> no 1's                         m1[nn, 2] := Vector(numelems(m1[nm2, 2])+1, {1=1});         # [1, 0, 0, ...]                         # [n-1,1] has 2 perms; adding 1's leaves nperms unchanged                                      nperms[nn, 2] := Vector([2, nperms[nm2, 2]]);                         # [n-1,1] has prod n-1; adding 1's gives prodnew = (prod + sum + 1) of p(n-2,2)                         prods[nn, 2] := Vector([n-1, prods[nm2, 2] +~ (n-1)]);                                         # k = 3;                         # appending 1 increases m1 by 1; adding 1's -> no 1's                         m1[nn, 3] := Vector([ m1[nm1, 2] +~ 1, Vector(numelems(m1[nm3, 3])) ]);                                 # appending 1: new nperms = (old nperms)*k/(m1+1); adding 1's leaves nperms unchanged                                                   nperms[nn, 3] := Vector([ zip((x,y) -> x*3/(y+1), nperms[nm1, 2], m1[nm1, 2]), nperms[nm3, 3]]);                         # [a,b] -> [a,b,1], p2 = (a*b) + (a+b); [a+1,b+1,c+1] has p2 = (a*b + b*c + a*c) + 2*(a+b+c) + 3                         p2[nn] := Vector([ prods[nm1, 2] +~ (n-1), p2[nm3] +~ (2*n-3) ]);                         # [a,b,1] has same prod as [a,b]; adding 1's gives prodnew = prod + p2 + sum + 1                         prods[nn, 3] := Vector([ prods[nm1, 2], prods[nm3, 3] +~ p2[nm3] +~ (n-2) ]);                 # calculate and store sum                 s := s + add(nperms[nn, 3] /~ prods[nn, 3]);                 ss[n] := s;                 # reindex Array rows for next iteration                 key := ListTools:-Rotate(key, 1);         end do;         ss end proc:
 > L2:=CodeTools:-Usage(sums2(200),iterations=4):

memory used=85.29MiB, alloc change=12.27MiB, cpu time=78.25ms, real time=452.50ms, gc time=3.91ms

 > EqualEntries(L,L2);

 >
 >

## solutions and export...

@raj2018 There are two pieces becouse there are two solutions for each kc and deltah, e.g., for kc=2.3, deltah = 0.8, M can be 0.0007832833567 or 1.134283416. If you don't want to see the bottom piece on the plot, just adjust the M range accordingly.

To get data out, it is best just to solve the equation for the values you want. An example is given in the attached

implicitplot3d.mw

## FormalPowerSeries...

@sursumCorda If you know the answer, you can cheat. Both convert(z*hypergeom([1/4, 1/2], [5/4], z**4), FormalPOwerSeries) and convert(EllipticF(z, I), FormalPowerSeries) give the same power series, verifying they are the same.

They are converted to apparently different expressions with integrals with convert(..., Int), so interconverting might be possible by rearranging, but it wasn't obvious to me.

## nice...

@acer Nice solution! I'm a bit surprised the Iterator method is so much slower than @sursumCorda's more-or-less equivalent method with combinat routines.

## eqns and variables...

Just to add to @mmcdara's comments. I especially agree with his last comment. Please clarify that you want to find each of  the 12 variables n1,n,2,n3,n4,4,a1,a2,a3,a4,s1,s2,s3,s4 in terms of the 4 parameters b,q,r,theta. So as noted you need 12 equations for your 12 variables.

The solutions for b=0 have arbitrary values of n3,n4 and s4. Do you expect this to be true in the general case? A dimensional analysis might help you to reduce the number of variables or parameters.

The reason it works for b=0 is that it can be cast as a polynomial system for that case. That can be done also for the full system, which might make it feasible to solve. But it would be better to simplify first, and we definitely need number of eqns = number of variables. Are the ai, ni and si components of 4-vectors? If so, perhaps the original equation(s) in matrix vector form would be useful.

## Events...

@dharr And here's how to do it with an event.

 > restart;
 > k12:=0.0452821; k21:=0.0641682; k1x:=0.00426118; Dose:=484; #Inj:=ifelse(t=1,Dose,0); ode_sys:=diff(SA2(t),t)=SA1(t)*k12-SA2(t)*k21,          diff(SA1(t),t)=-k1x*SA1(t); ics:=SA2(0)=0,SA1(0)=0;

 > Solution1:=dsolve({ode_sys,ics},{SA1(t),SA2(t)},type=numeric,events=[[t=1,SA1(t)=SA1(t)+Dose]]);

 > plots:-odeplot(Solution1,[[t,SA1(t),color="red",thickness=1],[t,SA2(t),color="blue",thickness=1]],t=0..50,labels=["Time t","Salicylic acid"],labeldirections=[horizontal,vertical]);

 >

## sort...

@vv Thanks, I just assumed that was the default order. As you undoubtedly know, the exact order isn't important, just that it is consistent for the left and right eigenvectors.

## yes...

@Christopher2222 Yes, you are right. Sorry for the confusion. The toggle input/output display works for plots and also if the output is on a different line because it was generated with "enter". But for equations where the input and output are on the same line, it seems you have to use view show/hidecontents and check/uncheck input or output.

## select input...

@Christopher2222 I only have 2024 right now but I think it is the same. After "plot(x^2)" then "ctrl=" you should have the plot on the same line, to the right of the command. Put the cursor in the plot command so that you are in the same document block (selecting the plot itself should also work). Then Edit->Document blocks->toggle input/output display makes the command disappear.

## bug...

@mmcdara local gamma; evalf(gamma); returns gamma and not the number in Maple 2024, i.e., the bug is fixed.

## simplify...

In general, use simplify(A-B) to ckeck if A and B are the same. But your worksheet has several problems. You cannot use N and N[1] as separate variables; once you assign to N[1], N no longer has the value you set previously. You can avoid this issue by using different names, which can be done by making subscripts by, for example, N__1. The same problem occurs with x and x[01], y and y[01] and maybe other places.

Do not use Digits:=1, which means to calculate with one digit accuracy. If you just want to display 1 digit, use Tools->Options->precision->round screen display.

You should not compare two indefinite integrals, since the arbritary constant could be different in the two cases.

@MaPal93 I didn't go through it all - just two places you asked questions about.

derivatives_final.mw

Two second derivatives change sign. For each, I find the zero with solve() but it returns two zeros instead of just one as I expected. Why?

I don't know - perhaps it is a solution on another branch of Lambda. Use fsolve if you just want a numerical answer.

I find minima and maxima with [DirectSearch]GlobalOptima, which is most of the time accurate but quite slow. Any alternative that is more efficient? I write "most of the time" because occasionally it gives me some very random results, completely off from the peaks and troughs I can eyeball from the plots.

GlobalOptima is overkill here - you only want a local maximum - Search is fast here, but needs an initialpoint. But here I would just use fsolve on the derivative with a range or initial guess.

My axes are always gamma*sigma_d^2 vs gamma*sigma_d*sigma_v. Are there any other param combinations that I can use? Would it make sense to scale the plots together with the actual axes, e.g., dividing both by gamma*sigma_d so that my new y-axis and x-axis are, respectively, sigma_d and sigma_v? How to do so? By doing so perhaps the plots are easier to interpret...

Yes, I wondered about that. Some people would non-dimensionalize the problem right at the beginning and then everything would be in terms of diminsionless quantities and the reader would be left to convert back to the real quantities. Back at the beginning you had f1 =Lambda*(sigma__v/sigma__d), so let's write a = sigma__v/sigma__d. and you also have Gamma = gamma*sigma__v*sigma__d, let's write b = gamma*sigma__v*sigma__d. So all quantities you write have to just have combinations of a and b. Here your axes are b/a and b. So for a universal plot with  axes sigma__d and sigma__v, those quantities would have to be expressible in terms of a and b, which I don't think is possible (you can't get rid of gamma unless you just take combinations of a alone, which can't be separated into sigma__v and sigma__d). Another way to think about this is to suppose you had such a plot, and you chose a point on the line for certain values of sigma__d and sigma__v. Does that result apply for all gamma? So the value of gamma seems essential to the problem, which may make sense to you in terms on the original problem.

## mathematical...

@sand15 I usually use it in combinatorial applications where being zero for b>0 makes mathematical sense. Wikipedia gives a "definition" that it is the coefficient in the expansion (1 + x)^alpha = sum(binomial(alpha,k)*x^k,k=0..infinity), for any complex alpha and |x|<1. Then  (1+x)^2 = 1*1 + 2*x + 1*x^2 + 0*x^3 + ...

DLMF defines it for complex alpha (=z) by 1.2.6, so binomial(2,3) = 2(2-1)..(2-3+1)/3! = 2(2-1)(2-2)/3! = 0

You are right about the conversions, but I think it is usual when converting between special functions that some cases might be exceptions.