pagan

5147 Reputation

23 Badges

17 years, 126 days

 

 

"A map that tried to pin down a sheep trail was just credible,

 but it was an optimistic map that tried to fix a the path made by the wind,

 or a path made across the grass by the shadow of flying birds."

                                                                 - _A Walk through H_, Peter Greenaway

 

MaplePrimes Activity


These are answers submitted by pagan

restart:

F:=proc(x::posint)
  local count, X;
  if x=1 then return [1,0]; end if;
  X := x;
  count := 1;
  while X<>1 do
    count := count+1;
    if X=2 then X:=1;
    elif irem(X,2)=0 then X:=X/2;
    else X:=3*X+1;
    end if;
  end do;
  [x,count-1];
end proc:

F(24);

(1)

st:=time():
plots:-pointplot([seq(F(i),i=1..10000)]);
time()-st;

 

(2)

g:=proc(x::posint,c)
  option remember;
  if x=1 then [1,c]; end if;
  if x=2 then [2,c+1];
  elif irem(x,2)=0 then g(x/2,c+1);
  else g(3*x+1,c+1);
  end if;
end proc:
G := proc(x)
  [x,g(x,0)[2]];
end proc:

G(24,0);

(3)

st:=time():
plots:-pointplot([seq(G(i),i=1..10000)]);
time()-st;

 

(4)

Download recurse.mw

Change w*-221.862 inside the assignment to sys to be -221.862*w instead. You need at least that much of a change, to correct the syntax error.

 

It wasn't quite clear to me whether you wanted to filter out a range of frequencies (possibly all above some value) or whether you wanted to allow only a range or band of frequences to get passed through. If you want to remove all frequences above a given value then either of those apporaches should suffice.

The examples below attempt to show a band getting passed through.

restart:


freqpass:=proc(M::Vector,rng::range(nonnegint))
local Mt, Xlow, Xhi;
  Xlow,Xhi := op(rng);
  Mt := DiscreteTransforms:-FourierTransform(M):
    if Xlow>1 then ArrayTools:-Fill(Xlow-1,0.0,Mt,0,1); end if;
    if Xhi<op(1,Mt) then ArrayTools:-Fill(op(1,Mt)-Xhi-1,0.0,Mt,Xhi+1,1); end if;
    Mt[1]:=Mt[1]/2;
    2*map(Re,DiscreteTransforms:-InverseFourierTransform(Mt));
end proc:

Make a sample signal M

f := x -> ((x-1/3)^5) + 50*sin(68*x):

a,b,n := -Pi,Pi,1000:
M := Vector(n,(i)->(evalhf@f)(a+i*(b-a)/n),datatype=float[8]):

plots:-pointplot([seq([k,M[k]],k=1..op(1,M))],style=line,
                  view=[0..1000,-600..300]);

 

filtered := freqpass(M,1..67):

plots:-pointplot([seq([k,filtered[k]],k=1..op(1,filtered))],style=line,
                 view=[0..1000,-600..300]);

 

filtered := freqpass(M,68..71):

plots:-pointplot([seq([k,filtered[k]],k=1..op(1,filtered))],style=line,
                 view=[0..1000,-600..300]);

 

filtered := freqpass(M,72..100):

plots:-pointplot([seq([k,filtered[k]],k=1..op(1,filtered))],style=line,
                 view=[0..1000,-600..300]);

 

f := x -> sin(5*x) + sin(118*x) - sin(130*x) + sin(123*x) - sin(127*x)
              - 1/3*sin(113*x) + 1/5*sin(300*x):

a,b,n := -Pi,Pi,1000:
M := Vector(n,(i)->(evalhf@f)(a+i*(b-a)/n),datatype=float[8]):

plots:-pointplot([seq([k,M[k]],k=1..op(1,M))],style=line,view=[0..1000,-7..7]);

 

filtered := freqpass(M,1..120):

plots:-pointplot([seq([k,filtered[k]],k=1..op(1,filtered))],style=line,view=[0..1000,-7..7]);

 

filtered := freqpass(M,1..100):

plots:-pointplot([seq([k,filtered[k]],k=1..op(1,filtered))],style=line,view=[0..1000,-7..7]);

 

 

Download fftfun00.mw

 

All this is doing is zeroing the result of the FFT except for the range you specify, and then doing the inverse-FFT. With a few tricks to try and stay efficient. Remember, the result of the FFT is (real part) frequency and (imaginary part) shift.

I'm not doing anything to allow you to express things in terms of time, hertz, sample duration. That's just bookkeeping and nomenclature. Bells and whistles. But let me know if this aspect has you puzzled. These samples below are just a number of data points, so how you interpret them according to some sampling speed, etc, is up to you. With this below, there are simply a number of discrete frequencies (depending on the number of data points and whatever you say the elapsed duration is) and you can pass them through or filter them out.

Notice the magnitude of the (supposedly non-noisy) lower frequency filtered result in the latter example. In this case the higher frequencies (now filtered out) contained a measurable portion of the signal strength.

If you do want to use such terminology then you have to say what you want done in the case that either frequency bound you specify doesn't match the granularity of your data. For given data and a given duration in units of time there will be a minimal attainable frequency determined by the number of nodes of data. If you specify a frequency cut-off that doesn't fall on a boundary then you have to agree how to handle it: eg. apportion it between the pair of neighbouring attainable discrete frequencies, or apportion it all to the single nearest attainable frequency.

You can issue unassign('t') or even t:='t' in order to undo the assignment to t.

But there is a way to manage your workflow which does not involve having to assign, unassign, reassign, and recompute with variables like t every time you want to make a new formula. Simply never assign to t at all! Just evaluate your formula at any new value of t when needed.

> formula := t^2+sin(s*t)+s^3;
                        2               3
                       t  + sin(s t) + s 

> eval(formula, t=2);
                                       3
                       4 + sin(2 s) + s 

> eval(formula, [t=2, s=3]);
                          31 + sin(6)

> list_of_values := [t=2, s=3]:

> eval(formula, list_of_values);
                          31 + sin(6)

> set_of_values := {t=3, s=7}:

> eval(formula, set_of_values);
                         352 + sin(21)

> new_formula := 11*t - cos(t);
                         11 t - cos(t)

Yet another approach is to create an operator or procedure from your formula, which can then be called with the desired value.

> restart:

> formula := t^2+sin(s*t)+s^3;
                        2               3
                       t  + sin(s t) + s 

> func := unapply(formula,[t,s]);
                            2               3
                 (t, s) -> t  + sin(s t) + s 

> func(2,s);
                                       3
                       4 + sin(2 s) + s 

> func(2,w);
                                       3
                       4 + sin(2 w) + w 

> func(2,3);
                          31 + sin(6)

This seems to be a FAQ. I guess that the manual and some tutorials encourage the assignment method, and then don't adequately explain why you might not want to do that.

You could try to contact him via his Mapleprimes profile. Using the contact link here should send an email to the address he has associated in his profile.

I think that those "oscillations" are artefacts of numeric float computation only, ie. noise.

Ok, so discont doesn't finish quickly/well when applied to diffTermScores[2] which is a shame. And fdiscont returns a whack of points. But we can determine just by looking at it that there are only a few kinds of piecewise in it. And they involve just a few key points on t=0..1.

> indets(diffTermScores[2],piecewise);

 /         /    1         1             1      1\  
{ piecewise|t < -, 0, t = -, undefined, - < t, -|, 
 \         \    2         2             2      3/  

           /    1  1      1             1       \           /
  piecewise|t < -, -, t = -, undefined, - < t, 0|, piecewise|
           \    2  3      2             2       /           \

      1         1                 4  1      4             4       
  t < -, 0, t = -, undefined, t < -, -, t = -, undefined, - < t, 0
      2         2                 5  3      5             5       

  \           /    1  -1      1                 4         4  
  |, piecewise|t < -, --, t = -, undefined, t < -, 0, t = -, 
  /           \    2  3       2                 5         5  

             4      1\\ 
  undefined, - < t, -| }
             5      3// 

So it's not hard to see where to break 0..1 apart so as to do a few separate `simplify` calls.

> newdTS2 := piecewise(
>    t<1/2,`assuming`([simplify(diffTermScores[2])],[t<1/2]),
>    t>1/2 and t<4/5,`assuming`([simplify(diffTermScores[2])],[t>1/2,t<4/5]),
>    t>4/5,`assuming`([simplify(diffTermScores[2])],[t>4/5]),
>    undefined);

         /    1       30        1             4       30        
piecewise|t < -, -------------, - < t and t < -, -------------, 
         |    2              2  2             5              2  
         \       (-11 + 10 t)                    (-11 + 10 t)   

  4                  \
  - < t, 0, undefined|
  5                  |
                     /

Great. You can plot or integrate that.

Now, presumably you want to figure out those problem points automatically. (I didn't look to see whether they even change with your Matrix, but if they don't then you're done now.)

This is a bit wordy, but seems to work

> convert(select(Y->type(Y,numeric),map(op,map(Y->map2(rhs@op,1,PiecewiseTools:-ToList(Y)),
>         indets(diffTermScores[2],piecewise)))),list);
                             [1  4]
                             [-, -]
                             [2  5]

Suppose you had a longer list of key points. You could work on it like this.

> U:=[1/3,1/2,2/3,4/5];
                          [1  1  2  4]
                          [-, -, -, -]
                          [3  2  3  5]

> t<U[1],G,seq(op([t>U[i-1] and t<U[i], G]), i=2..nops(U));
           1     1             1     1             2     
       t < -, G, - < t and t < -, G, - < t and t < -, G, 
           3     3             2     2             3     

         2             4   
         - < t and t < -, G
         3             5   

Which leads me to this.

> U:=convert(select(Y->type(Y,numeric),map(op,map(Y->map2(rhs@op,1,PiecewiseTools:-ToList(Y)),
>                   indets(diffTermScores[2],piecewise)))),list);

                             [1  4]
                             [-, -]
                             [2  5]

> piecewise(t<U[1],`assuming`([simplify(diffTermScores[2])],[t<U[1]]),
>           seq(op([t>U[i-1] and t<U[i],
>               `assuming`([simplify(diffTermScores[2])],[t>U[i-1] and t<U[i]])]),
>               i=2..nops(U)),undefined);

         /    1       30        1             4       30        
piecewise|t < -, -------------, - < t and t < -, -------------, 
         |    2              2  2             5              2  
         \       (-11 + 10 t)                    (-11 + 10 t)   

           \
  undefined|
           |
           /

Now someone will show a much easier way. ;) What happens if you Matrix data is exact rational and not floats to begin with?

The problem is that some of the subexpressions in your big expression T are actually atomic identifiers.

An atomic identifier is a Maple name, formed of single left-quotes between which lies magic stuff. Maple's Standard GUI knows how to prettyprint the stuff as if it were some 2D Math. But underneath the whole subobject (whether it be a complicated fraction or formula or whatever) is just a single name.

I'd like to know how those subexpressions became atomic identifiers. If it happened automatically, without a deliberate action by you, then it may be a bug (because the imaginaryunit was set?).

For example, your expression T has what looks like the following subexpression in it.

R/(j*omega*C[s])

But that subexpression is not what it seems. In fact, the object T contains just a fancy name which gets type-set as that. Underneath, the subexpression is actually the name

        `#mi("R")`/(j*omega*`#msub(mi("C"),mi("s"))`)

As someone else mentioned, you can try using the right-click context-menus to convert all of T at once from 2D Math to 1D Maple notation. That is, select the input T (not the output), right-click and make the context-menu choice 2D Math -> Convert To -> 1D Math Input. That's probably easiest, given the fact that you've already entered it in 2D Math.

Alternatively, you could try and replace each of the atomic identifier subexpressions inside T, using repeated application of set of similar context-menu actions. But to do so with the mouse you'd first have to find which subexpressions in particular are atomic identifiers. That is to say, how would you know how much to select with the mouse, before right-clicking? You could lprint(T), and try and interpret it visually, but it's a little messy. Or you could try and remove all atomic identifiers with a program.

> FromAtomic := proc(expr)
> local temp, origkom;
>   origkom := kernelopts(opaquemodules=false);
>   try
>     if assigned(Typesetting:-SemanticTable[convert(expr,string)]) then
>       temp := Typesetting:-SemanticTable[convert(expr,string)];
>       temp:=op(-1,Typesetting:-Parse(temp));
>     end if;
>   catch:
>   finally
>     kernelopts(opaquemodules=origkom);
>   end try;
>   if assigned(temp) then
>     eval(temp);
>   else
>     expr;
>   end if;
> end proc:

> subsindets(T,name,FromAtomic);
           //         /         I    \ /
            |omega Cp |R L - --------| |
  - (I R L) |         \      omega Cp/ |
            |                          |
            \                          \
                I R                         I R L          \\
  - --------------------------- - -------------------------||
               /        I     \            /         I    \||
    omega C[s] |R - ----------|   omega Cp |R L - --------|||
               \    omega C[s]/            \      omega Cp///

> simplify(%);

(-R omega C[s] + I) L - --------------------------------------- R L omega Cp - I + L R omega C[s] - I L

Anyone interested in "robustifying" that `FromAtomic` routine might have a peek at Typesetting:-RemoveAtomic. But you can't just utilize that routine as is, since it has the here unwanted side-effect of not just returning the non-atomic expression but also of undesirably unassigning the associated entry from the Typesetting:-SemanticTable.

Let's leave aside programming methodology for now. (There are some easier ways to construct your procs and do mapped differentiation.)

Try this change:

fDiffx:=proc(X,methodList)  # differentiate functions based on t
	local i,tmp;
	for i in methodList do:
           tmp[i]:=diff(simplify(convert(X[i],rational)),t) assuming t>0, t<1;
	od;
	Vector(3,i->tmp[i]);
end proc:

And this one (so that it does only numeric integration, and doesn't "go to town" trying any symbolic integration):

distortion:=evalf(Int(tmp,0.0..1,epsilon=1e-4,digits=10));

I got a result of 0.2453926822 for that last one.

plot(tmp,0..1);

Without the convert,rational in fDiffx (but with the `simplify` and `assuming`) I had a problem with the first part of that sawtooth graph, where it computed with lots of *numerical* discontinuties.

Of course, I also changed a line as

diffTermScores:=fDiffx(termScores,[1,2,3]):

You want to eliminate intermediate variables? One simple approach could involve calling the `eliminate` command. (How you would further analyze its output for more complicated examples is more interesting.)

restart:
sys:={y = z1 + z2^2,
      z1 = z2 + 3*exp(u),
      z2 = 2*u^3*z1}:
T:=eliminate(sys,[z1,z2]):
new:=normal(isolate(op(T[2]),y)):
expected:=y=3*exp(u)/(1-2*u^3) * (1 + 4*u^6*(3*exp(u)/(1-2*u^3))):
evalb(normal(rationalize(new))=normal(rationalize(expected)));
                              true

You could contact Maplesoft's Tech Support. But for mild excitement, maybe you could first try using chroot for running one of the instances.

You can do either of these

for i from 1 by 1 to 5 do Y[i] :=i ; X[i] :=i^2 end do:

for i from 1 by 1 to 5 do (Y[i],X[i]) :=i, i^2 end do:

Some of your expressions are formed with the original assumption on u1. But then you place new assumptions on u1 when you do the three lines starting with assume(uSTARm2>=...). And then you form new expressions using that u1 which has now got different assumptions. Now you have two variables u1, each with different assumptions on them. It would be the same even if you used `additionally` in the second group. This is a common muddle that people get into, with assumptions.

Why not just remove ALL the lines which call `assume`, and then do

   fsolve(simplify(e4)) assuming u1>0;

Sometimes you can use the Array constructor to make such tasks convenient. So perhaps that can help you with this (or a future similar) task where you have a variable number of arguments.

Also, isn't your problem symmetric in all the q[i]? I mean, s[3,2,1,3,2] will be the same as s[1,2,2,3,3], yes? If so then you can use a symmetric indexing function to avoid repeating work, although at the cost of sorting indices in the implementation below.

> restart:

> `index/commutingindexes`:=proc(idx,M,val)
> local idx1:
>   idx1:=op(sort(idx));
>   if nargs = 2 then M[idx1];
>   else M[idx1]:=op(val); end if:
> end proc:

> makeT := n -> Array(commutingindexes,seq(1..3,i=1..n),()->`+`(args)):

> Z := makeT(5):
> Z[1,2,2,3,3];
                               11
> Z[2,1,3,2,3];
                               11

> makeS := proc(n)
> local S, i, summands, k, j;
>   S:=table():
>   for i from 0 to 3^n-1 do
>      summands:=NULL:
>      k:=i:
>      for j to n do       
>         summands:=summands, (k mod 3) +1:
>         k:= floor(k/3):
>      end do:
>      S[summands]:=convert([summands],`+`):
>   end do:
>   S;
> end proc:

> seq( print(time(makeS(ii))/time(makeT(ii))), ii=8..12 );
                          14.04678363
                          12.59624413
                          17.11492063
                          17.12746191
                          17.53507271

I'm not suggesting using an Array as a best approach. But if the number of variables (dimension) is going to become very high then perhaps any savings is good.

Another question: are you planning on using all entries of the table, or most, or just a few? If just a few, then could you not just implement a function to do it on demand? Eg,

> H:=proc() `+`(args); end proc:

> H(3,2,1,1);
                               7

A question for the experts: is there any easy way to get a symmetric remember table? I ask because unlike rtables the 'symmetric' indexing function of "older" tables works with more than 2 dimensions (and might be more efficient that the crude `sort` inside the `index/commutingindices` routine at top).

> H := subsop(4=table([],symmetric),eval(H)):
Error, remember table can only be set to NULL

Are you still using Classic or the Standard GUI, in 14? It wasn't clear from your post.

Perhaps you might even have thought that 14 had no Classic, if you didn't get a DeskTop icon for it when you installed. If that's the case, and you have 32bit Maple 14 installed, then look for in under the Start Menu or in the bin.win Maple 14 folder as the cwmaple.exe file.

Or perhaps you are already using Classic in 14 from 32bit Maple for Windows, in which case you should send the problematic worksheets to Maplesoft's Tech Support.

Or perhaps you are using Windows XP/Vista 64, and 64bit Maple 14 for which there is no distributed Classic for Windows.

Or, if you are using Standard (intentionally, and intend on continuing to give it a shot) then have you tried saving your worksheets as .mw instead of as .mws? (I mean, maybe the issue is with Standard trying to auto/save as old .mws format.)

Using with(MyLib) as a way to resolve the name `h` as used inside the already authored or loaded package AnotherPackage is bad form (and this is an example of why). Also, it really matters which version (global, or export of MyLib) of the name `h` is bound when AnotherPackage is defined and saved,

I suggest that you write AnotherPackage using the syntax MyLib:-h() inside the code of AnotherPackage, whenever you need to call export `h` of MyLib.

If done that way, properly, you shouldn't even have to with(MyLib) in order to `with` or use AnotherPackage. You'd just have to make sure that MyLib is accessible in the libname path whenever you used AnotherPackage.

You might also utilize `use`, but there are some caveats abour which form works I think (and `uses` might not be ok at all).

First 27 28 29 30 31 32 33 Last Page 29 of 48