Carl Love

Carl Love

28095 Reputation

25 Badges

13 years, 100 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

If you have any procedure, and you want to expand its domain to types that it wasn't designed for, or you want to change its behavior for certain types of input, then you can use overload. In this case, we want to extend the definition of numelems to include input of type name.

unprotect(numelems);
numelems:= overload([
     proc(C::name &under (C-> eval(C))) option overload; 0 end proc,
     numelems
]):

numelems(a);
     0

numelems([1,2,3]);

     3

numelems(a+b);

Error, invalid input: numelems expects its 1st argument, t, to be of type indexable, but received a+b

(A curiosity about the above code: I was forced to use C-> eval(C) rather than the simpler eval. The latter tries to use two-argument eval, as seen in the following:

type(a, name &under eval);

Error, (in type/&under) invalid input: eval received NULL, which is not valid for its 2nd argument, eqns

This is a very curious error message because I thought that a procedure parameter could never be passed NULL! But what I'm more curious about is why it chooses to use two-argument eval. Here's a simpler example of the same thing:

eval(a, NULL);

Error, invalid input: eval received NULL, which is not valid for its 2nd argument, eqns)

Your entire loop needs to be in one execution group. To achieve this, use the editing command Join Execution Groups. You can find this by going to the Edit menu, then Split or Join, then Join Execution Groups. The shortcut key for this is F4. So, the entire operation can be done by placing your cursor at the top of the worksheet and repeatedly pressing F4. The cursor won't move, but if you look closely at the execution group boundaries at the left edge of the worksheet, you'll see them change.

After having merged the execution groups, add the loop header to the top and end do; to the bottom. If any of your original execution groups contained text rather than code, it doesn't matter: The parser will treat that text as comments. So, this process also gives you a way to include formatted text comments in your code.

This is the same error message as you got with your bubble sort procedure a few days ago. Like I told you, you can't ordinarily assign to a formal parameter, such as n. In other words, n shouldn't appear on the left side of :=. If you need to do something like that, first copy n to a local variable, as Markiyan has copied it to k. Also, I told you not to use the print statement as a substitute for a procedure's return value.

Here's my Collatz procedure:

Collatz:= proc(N::posint)
local n:= N, count;
     for count while n <> 1 do
          n:= `if`(n::even, iquo(n,2), 3*n+1)
     end do;
     count
end proc:

I find binary infix operators very useful. In this case, I agree with the use of the notation of that "other CAS", although I usually abhor it. The following code allows for the use of notation very similar to what you posted. You just need to change the angle symbol to &<.

phasor:= module()
option package;
export
     `+`:= proc(Z1::`&<`(realcons,realcons), Z::seq(`&<`(realcons,realcons)))
     option overload;
     local z, r:= add(op(1,z)*exp(op(2,z)*Pi/180*I), z= args);
          evalf(abs(r)) &< evalf(argument(r)*180/Pi)
     end proc,

     `-`:= proc(Z::`&<`(realcons,realcons))
     option overload;
          (-1*op(1,Z)) &< op(2,Z)
     end proc
;
end module:

with(phasor):

4&<45 + 5&<30;

This allows for sums and differences with an arbitrary number of phasors. With a little more work, I could get the output to print without the & and the extra parentheses.

After you fix the error pointed out by Robert Lopez, you should get an error about pi. In Maple, this must be capitalized (Pi) if you mean the famous mathematical constant.

Here's an example of what you want. It involves using the keyword unevaluated function typeset in the plot command (see ?plot,typesetting). I tested with 99 frames of triplets of random numbers, and there's no jitter. All that you need from this are the procedures P3 and T3; the rest is just code I wrote to test them.

P3:= (x::realcons)-> sprintf("= %5.3f   ", x):
frames:= 99:
L||(1..3):= 'RandomTools:-Generate(list(float(method= uniform), frames))' $ 3:
T3:= (f::seq(realcons))-> typeset(seq([theta[_k], P3(f[_k])][], _k= 1..nargs)):
plots:-display(
     [seq(
          plot(
               0,
               title= T3(seq(L||_k[_j], _k= 1..3)),
               titlefont= [TIMES,ROMAN,14]
          ), _j= 1..frames
     )],
     insequence
);


It can be done with a single command:

a2_||(pos,neg):= selectremove(t-> sign(t) = 1, a2);

Clearly you want a matrix whose entries are matrices rather than a block matrix because your sizes won't fit into a block matrix. Here's how to do it:

a1 := Matrix(3, [1, 2, 3, 7, 8, 9, 13, 14, 15]);
a2 := Matrix(3, 2, [5, 6, 11, 12, 17, 18]);
a3 := Matrix(2, [19, 20, 25, 26]);
a4 := Matrix(3, 2, [5, 6, 11, 12, 17, 18]);
A := Matrix(2, 2, (i,j)-> a||(2*(i-1)+j) );

Here's a complete example showing the method that I outlined in my Replies. I compute, by purely numeric computation, the surface area of the isosurface MF(x,y,z) = 11 for MF(x,y,z) = sin(y*z)+x+11. Then, for accuracy comparison, I compute the surface area by trivial numero-symbolic computation.

My worksheet could be made more efficient; there's some redundant calculation. I used the way that I did for clarity of exposition.

 

restart:

Digits:= 15:  #Usually the best value for numerics.


(* The MF below is kind of trivial, but it's a numeric black box to
simulate what you have. It only returns numeric information---none of
the root-finding, derivative, or integration routines can see inside it
to analyze the symbolics. *)

MF:= proc(x,y,z)
option remember;
     `if`(andmap(type, [args], realcons), sin(y*z)+x+11, 'procname'(args))
end proc:

(* F is the procedure that expresses the isosurface MF(x,y,z) = 11 as a
function x = f(y,z) *)

F:= proc(y,z) option remember; fsolve(x-> MF(x,y,z)-11, -1..1) end proc:

(* Fy and Fz are procedures for the numeric computation of the partial
derivatives of F. *)

Fy:= proc(y,z)
local x;
     if andmap(type, [args], realcons) then
          x:= F(y,z);
          -fdiff(MF, [2], [x,y,z])/fdiff(MF, [1], [x,y,z])
     else
          'procname'(args)
     end if
end proc:

Fz:= subs([2]= [3], eval(Fy)):
 


#Now do a pure numeric computation of the surface area of the isosurface.
evalf(Int((y,z)-> sqrt(1+Fy(y,z)^2+Fz(y,z)^2), [0..1, 0..1], epsilon= 1e-7));

1.24007661357377

(1)


(* For accuracy comparison, use numero-symbolic computation to get
the same surface area. Simply use the method from elementary
multivariable calculus. *)
G:= (y,z)-> sin(y*z):

evalf(Int(sqrt(1+D[1](G)(y,z)^2+D[2](G)(y,z)^2), [y= 0..1, z= 0..1]));

1.24007661358042

(2)

#So the first integration had 11 digits accuracy even though we only asked for 7!
#Amazing!!

``

 

Download numeric_isosurface_surface_area.mw

Your immediate problem is that X is a formal parameter, and you can't assign to a formal parameter. It may appear to you that you are only assigning to the elements of X, but since it's a list, any assignment to it reassigns the whole list. This makes assigning to the elements of a list a very bad thing to do, even though Maple will sometimes allow you to get away with it. The solution is to convert the list to a Vector on input, do your sorting, and then convert back to a list on output.

Some other points:

  1. Sorting the empty list should return the empty list; it shouldn't be an error. In the code below, the empty list just "falls through" because the outer for loop will just be skipped if n is 0 or 1.
  2. You can switch two elements in a Vector in a single assignment without using an intermediate variable. Look carefully at the central assignment in my code below.
  3. Don't use the print statement to return values. Either use a return statement, or just put the expression to be returned immediately before the end proc.

Here is your code with all of these recommendations implemented:

Bubble := proc(L::list(realcons))
local n := nops(L), X := Vector(L), i, j;
     for i to n-1 do
          for j to n-i do
               if X[j+1] < X[j] then
                    X[j..j+1]:= X[[j+1, j]]
               end if
          end do
     end do;
     convert(X, list)
end proc:

#Test it.
R:= rand(-99..99):
L:= ['R()'$99];

     output omitted

Bubble(L);

     output omitted

All in all, you did a good job! Your logic was perfect, and your only real problem was due to an idiosyncracy in the way that Maple implements lists.

 

 

I followed the link from the Wikipedia page to the paper that has the Matlab code. Then I meticulously translated that into efficient and compact Maple code. The totality of the Matlab code accounts for only the GRF:= ...; paragraph of my code! The result is a greyscale image. To add color, I put the high values on a blue scale and the low values on a red+green scale.

restart:

Digits:= 15:
GRF:= proc(n,r)
uses AT= ArrayTools, IT= ImageTools;
local
     GRF:= (IT:-FitIntensity@IT:-Convolution)(
          AT:-RandomArray(n, 'distribution'= 'normal'),
          Matrix(2*r+1$2, (i,j)-> `if`((i-r-1)^2 + (j-r-1)^2 <= r^2, 1, 0))        
     )
;     
     (IT:-Preview@IT:-CombineLayers)(
          IT:-FitIntensity(map[evalhf](x-> `if`(x >= .5, 0, x), GRF))$2,
          IT:-FitIntensity(map[evalhf](x-> `if`(x < .5, 0, x), GRF))
     )
end proc:

GRF(300, 13);

When the number of parameters of a BVP is greater than or equal to the number of conditions at one of the boundary points, it can be solved as an IVP with fsolve being used to determine the parameters. I like this idea because I think that the main IVP solver (rkf45) is more robust than any of the BVP solvers. The following worksheet produces a very high-accuracy solution in good time using little memory, although I haven't compared it timewise or memorywise to Preben's. Accuracywise, my value of omega computed at Digits = 10 has three more digits of accuracy than Preben's. (I don't mean to denigrate his solution, for he just used the default error control.)

restart:


ODE:= diff(y(x),x$2) = -(8*omega*(1-exp(-8*x))*exp(-8*x)/2/x^2-Pi^2/32/x)*y(x);
ICs:= y(0)=0, D(y)(0)=1:
BC:= y(1/2)=0:
IVP:= unapply({ODE,ICs}, omega):

diff(diff(y(x), x), x) = -(4*omega*(1-exp(-8*x))*exp(-8*x)/x^2-(1/32)*Pi^2/x)*y(x)

(1)


#Option 'params' isn't allowed for initially singular IVPs, so we need
#a proc.
F:= proc(omega)
     Digits:= Digits+2;
     if omega::numeric then
          dsolve(
               IVP(omega), numeric,
               output= Array([op(lhs(BC))]),
               relerr= 10^(2-Digits)
          )[2][1][1,2] -
          rhs(BC)
     else
          'procname'(omega)
     end if
end proc:   

Sol:= CodeTools:-Usage(fsolve(F(omega), {omega=0}));

memory used=87.39MiB, alloc change=59.60MiB, cpu time=780.00ms, real time=785.00ms

 

{omega = .8054765485}

(2)

Digits:= 15:

Sol:= CodeTools:-Usage(fsolve(F(omega), {omega=0}));

memory used=101.00MiB, alloc change=43.60MiB, cpu time=858.00ms, real time=854.00ms

 

{omega = .805476548534690}

(3)

SolIVP:= dsolve(IVP(eval(omega, Sol)), numeric, relerr= 10^(2-Digits)):

plots:-odeplot(SolIVP, [[x,y(x)],[x,D(y)(x)]], 0..1/2);

 

NULL

 

Download BVP_as_IVP.mw

 

If I understand you correctly, you want to increase the voltage, V, by steps and do fsolve at the increased voltages, at each step using as initial values the solutions from the previous step. Is that correct? In the worksheet below, I've done all this using a step of 0.1. The keys to making this work is to leave symbolic and to use unapply to turn f1 and f2 into procedures with parameter V. The other changes that I made are just to keep everything neat and symbolic until it needs to made numeric.

restart:


#I haven't included V in the Params.
(V0, Vstep, Vend):= (10.0, 0.1, 11):
Params1:= [
     a1= 3e-6, a2= 42e-6, a3= 50e-6, d= 2.75e-6, b= 1000e-6, Kt= 1.54e-9,
     Ky= 6.49, m= 4.3e-11, IIm= 2.5e-20, tetamax= d/a3, e= 8.954e-12
]:
Params2:= [
     etay= IIm*e*b*V^2/2/Kt/d^2/m/tetamax, etaa= e*b*V^2/2/tetamax^3/Kt,
     alpha= a1/a3, beta= a2/a3, wt= sqrt(Kt/IIm), wy= sqrt(Ky/m)
]:
Params3:= [w0= wy/wt]:


#The next two lines are for display purposes only.
op~([Params||(2..3)]);
<op~((eval['recurse']@op)~([[Params2, Params1], [Params3, op~([Params||(1..2)])]]))[]>;

[etay = (1/2)*IIm*e*b*V^2/(Kt*d^2*m*tetamax), etaa = (1/2)*e*b*V^2/(tetamax^3*Kt), alpha = a1/a3, beta = a2/a3, wt = sqrt(Kt/IIm), wy = sqrt(Ky/m), w0 = wy/wt]

Matrix(7, 1, {(1, 1) = etay = 0.406358968726833e-2*V^2, (2, 1) = etaa = 0.174734356552538e-1*V^2, (3, 1) = alpha = 0.600000000000000e-1, (4, 1) = beta = .840000000000000, (5, 1) = wt = 248193.4729, (6, 1) = wy = 388497.4035, (7, 1) = w0 = 1.56530064618808})

D1:= 1-x3-beta*x1:  D2:= 1-x3-alpha*x1:
f1:= -x1+etaa/x1^2*((1-x3)*(1/D1-1/D2)+ln(D1/D2));
f2:= -w0^2*x3+etay/x1*(1/D1-1/D2);

-x1+etaa*((1-x3)*(1/(-beta*x1-x3+1)-1/(-alpha*x1-x3+1))+ln((-beta*x1-x3+1)/(-alpha*x1-x3+1)))/x1^2

-w0^2*x3+etay*(1/(-beta*x1-x3+1)-1/(-alpha*x1-x3+1))/x1


#The next line is the key step: Converting f1 and f2 into procedures with
#parameter V.
(f1,f2):= unapply~(eval['recurse']([f1,f2], op~([Params||(1..3)])), V)[]:


#The main computation:
R:= fsolve({f1,f2}(V0), {x1,x3}):

for V from V0 by Vstep to Vend do
     Sol['V'=V]:= R;
     R:= fsolve({f1,f2}(V+Vstep), R)
end do:


#The following line is for display purposes only.
<zip(`::`,(([indices],[entries])(Sol, 'indexorder', 'nolist')))[]>;

Vector(11, {(1) = (V = 10.0)::{x1 = 1.70749422232618, x3 = .989389491442727}, (2) = (V = 10.1)::{x1 = 1.71501418386313, x3 = .990235617861761}, (3) = (V = 10.2)::{x1 = 1.72250003242707, x3 = .991087061384685}, (4) = (V = 10.3)::{x1 = 1.72995227916113, x3 = .991943727104430}, (5) = (V = 10.4)::{x1 = 1.73737142358238, x3 = .992805521938024}, (6) = (V = 10.5)::{x1 = 1.74475795392970, x3 = .993672354583526}, (7) = (V = 10.6)::{x1 = 1.75211234749862, x3 = .994544135478267}, (8) = (V = 10.7)::{x1 = 1.75943507096395, x3 = .995420776758373}, (9) = (V = 10.8)::{x1 = 1.76672658069047, x3 = .996302192219486}, (10) = (V = 10.9)::{x1 = 1.77398732303236, x3 = .997188297278666}, (11) = (V = 11.0)::{x1 = 1.78121773462185, x3 = .998079008937410}})

 

``

 

Download 2DOF_mod.mw

There are various uses for ~ in Maple. Your syntax is wrong. For elementwise operation, the ~ goes after the procedure name. So that should be int8~(round~(S)).

Simply take your first example, and change type~(r,anything) to evalb~(r).

First 248 249 250 251 252 253 254 Last Page 250 of 395