## 13776 Reputation

14 years, 59 days

## Something like...

the attached, where I have picked arbitrary values for the parameters 'b' and 'c', since you didn't supply any.

 > restart: # # Pick arbitrary values for b and c since OP # does not supply these #   pars:= [b=4, c=6];   f:= unapply(eval((2*z^3+b*z^2+c)/(3*z^2+2*b*z), pars), z);   H:= Array(0..100):   H[0]:= 222.07;   for j from 1 by 1 to 100 do       H[j]:= f(H[j-1]);   od:   H;
 (1)
 >

## Well...

the attached seems to work

 > restart;   expr:= sqrt(alpha^2+k^2)*(-sqrt(alpha^2+k^2)*exp(k*z0)+k*exp(sqrt(alpha^2+k^2)*z0))          /          ((-alpha^2-k^2)*exp(k*z0)+exp(sqrt(alpha^2+k^2)*z0)*sqrt(alpha^2+k^2)*k);   simplify(expand(expr));
 (1)
 >

## The plot() command...

has a size=[w,h] option.

## After fixing...

some (minor) syntax errors, the basic problem with your code is that the section

```  for i from 2 to np-1 do
h_n[i]:=(i-1)*Hc/(np-1):
H:=h_n[i]:
theta_minus:= fsolve(equ,theta=theta_s..Pi):
m_n[i]:=evalf(cos(theta_minus-psi)):
end do:
```

fails to produce numeric values for all m_n[i]. This happens because, for i>=59, the equation 'equ' does not have a root in the specified range theta_s..Pi. The nearest root to the specified range is slightly to the left of theta_s. I have no idea what you want to do for the case where no root occurs. As a (very crude) workaround, the attached computes all roots in the range 0..Pi, then selects the one  which is closets to Pi.

Although much slower than using fsolve() with a correct range, at least this executes and results in a plausible graph - see attached

 > restart:   with(plots):   with(plottools):   w:=-Ps/2*(Hk*cos(theta)^2+2*H*cos(theta-psi)):   dw:=diff(w,theta):   equ:=dw=0;   Hg:=(sin(psi)^(2/3)+cos(psi)^(2/3))^(-3/2):   psi:=35.0*Pi/180:   H:=0.1:   Ps:=1.0:   Hk:=1.0:   Hc:=evalf(Hk*Hg):   theta:='theta':   theta_s:=evalf(Pi-arctan(tan(psi)^1/3)):   H:='H':   np:=100:   H_max:=1.0:   h_p:=Array(1..np):   h_n:=Array(1..np):   m_p:=Array(1..np):   m_n:=Array(1..np):   h_p[1]:=0.0:   m_p[1]:=cos(psi):      for i from 2 to np do       h_p[i]:=(i-1)*H_max/(np-1):       H:=h_p[i]:       theta_plus:=fsolve(equ,theta=0..Pi/2):       m_p[i]:=evalf(cos(theta_plus-psi)):   end do:   h_n[1]:=0.0:   m_n[1]:=evalf(cos(Pi-psi)):   for i from 2 to np-1 do       h_n[i]:=(i-1)*Hc/(np-1):       H:=h_n[i]: # # No root in range theta_s..Pi for i>=59 # #      theta_minus:= fsolve(equ,theta=theta_s..Pi):       theta_minus:= max(Student[Calculus1]:-Roots(equ,theta=0..Pi));       m_n[i]:=evalf(cos(theta_minus-psi)):   end do:   h_n[np]:=Hc:   m_n[np]:=evalf(cos(theta_minus-psi)):   pl_pp:=Array(1..np):   pl_pn:=Array(1..np):   pl_np:=Array(1..np):   pl_nn:=Array(1..np):   i:='i':   for i from 1 to np do       pl_pp[i]:=point([h_p[i],  m_p[i]],color=red):       pl_np[i]:=point([h_n[i],  m_n[i]],color=blue):       pl_pn[i]:=point([-h_n[i],-m_n[i]],color=red):       pl_nn[i]:=point([-h_p[i],-m_p[i]],color=blue):   end do:   display( [ seq(pl_pp[i],i=1..50),              seq(pl_np[i],i=1..50),              seq(pl_pn[i],i=1..50),              seq(pl_nn[i],i=1..50)            ]          );
 >

## Some observations...

```  expr:=(s+exp(-Pi*s)-exp(-2*Pi*s))/(s*(s^2+2*s+2));
H:=inttrans:-invlaplace(expr,s,t);
Y:=convert(H,piecewise);```

using the "trick" _EnvUseHeavisideAsUnitStep := true, changes the definition of the function.

Instead of having (the correct) Heaviside(0)=undefined, so that Heaviside(t-Pi) and Heaviside(t-2*Pi) lead to undefined results for t=Pi and t=2*Pi, you are specifying that Heaviside(0)=1, so that the underlying function is now defined for all t - that is a different function!!

` f := piecewise(v < 0, -c, v=0, 0, v > 0, c);`

the function 'f' is defined for all 'v'. Converting to a Heaviside representation means that you will be introducing a value of 'v' for which the function is undefined - so again you are changing the function definition. Like the help for convert(expr, Heaviside) says (emphasis added)

Note that conversion to Heaviside can result in a function which differs from the original in finitely many points.

I don't think I have ever used _EnvUseHeavisideAsUnitStep := true - too risky!!

I have occasionally used the NumericEventHandler() as in the following -

```  restart;
f := piecewise(v < 0, -c, v=0, 0, v > 0, c);
NumericEventHandler(invalid_operation = `Heaviside/EventHandler`(value_at_zero = 1/2)):
convert(f, Heaviside);
eval(%, v=0);
```

but I pretty much always set it so that Heaviside(0)=1/2 - which just somehow seems a bit more "symmetrical", (and for the above example, gives the correct result).

## More roots...

can be obtained, as shown in the attached. Note that the function is even, so for every positive root there is a negative root.

 > restart:   Digits:=20:   expr:=cos(x)*cosh(x)-1;   plot( expr, x=0..20, view=[default, -10..10]);   rts:=Student[Calculus1]:-Roots(expr, x=0..20, numeric=true);
 (1)
 >

## It is not too difficult...

to come up with a simple function which will perform the required operation, on all the examples you provide. (see the attached.)

However, this "simple" function will start to get (a lot!) more complicated for example, if it has to recognise that

2*f(a)+f(b )= f(a)+f(a)+f(b ) --> f(2*a+b)

or indeed handle any term of which the target function is just a part, for example, would you expect

(2+sin(x))*f(a)+f(b) = 2*f(a)+sin(x)*f(a)+f(b)=f(a)+f(a)+sin(x)*f(a)+f(b) --> sin(x)*f(a)+f(2*a+b)

or something else??

Maybe you can guarantee that these more general cases never occur?

 > restart;   expr1:= f(A) + f(B);   expr2:= f(A) + f(B) + f(c);   expr3:= f(A) + f(B) + f(c) + f(10*c);   expr4:= sin(x) + f(A) + f(B) + 10*exp(x)/13 + cos(x) + f(c) + f(10*c);   func:= (z, ee)-> local j; z( add                                ( op(1, j),                                  j in select                                       ( j->evalb( op(0,j)=z),                                         ee                                       )                                )                              )                              +                              remove                              ( j->evalb( op(0,j)=f),                                ee                              ):   func( f, expr1);   func( f, expr2);   func( f, expr3);   func( f, expr4);
 (1)
 >

## Use a dsolve(....,numeric)...

as in the attached.

 > restart;   with(plots):   ode:= diff(y(x), x\$5)=exp(-x)*y(x)^3;   bcs:= y(0)=1, D[1](y)(0)=1/2, D[1,1](y)(0)=1/4, y(1)=exp(1/2), D[1](y)(1)=exp(1/2)/2;   sol:= dsolve( [ode, bcs],                 numeric               ):   odeplot( sol,            [x,y(x)],            x=0..1          );
 >

## Interesting one...

First of all there is no generally agreed definition of "autocorrelation" - see https://en.wikipedia.org/wiki/Autocorrelation

Maple has two (differently-defined) autocorrelation commands, namely

SignalProcessing:AutoCorrelation()
Statistics:-AutoCorrelation()

These two commands give somewhat different results, and you "appear" to be using the latter (see the attached). Before you proceed much further, I can only suggest that you read the help for both of these commands to work out which definition of autocorrelation is the closest for your needs (even then you may have to do a bit of 'scaling' etc).

 > restart;   V:=Vector([seq( evalf(sech(j-2048)), j=0..4096)]): # # Plot the signal #   plots:-listplot(V, view=[1850..2200, default], style=pointline, color=red); # # Perform autocorrelation from SignalProcessing() package #   AC1:=SignalProcessing:-AutoCorrelation( V ): # # check min/max values in autocorrelation #   min[index](AC1);   AC1[%];   max[index](AC1);   AC1[%]; # # Plot autocorrelation signal, and zoomed plot for small 'y-value' #   plots:-listplot(AC1, style=pointline, color=red);   plots:-listplot(AC1, view=[default, -0.01..0.01], style=pointline, color=red);
 > # # Perform autocorrelation from Statistics() package #   AC2:=Statistics:-AutoCorrelation( V , 0..4096): # # check min/max values in autocorrelation #   min[index](AC2);   AC2[%];   max[index](AC2);   AC2[%]; # # Plot autocorrelation signal, and zoomed plot for small 'y-value' #   plots:-listplot(AC2, style=pointline, color=red);   plots:-listplot(AC2, view=[default, -0.01..0.01], style=pointline, color=red);
 >

## Trivial...

See the help at ?Export, or the attached - where you will have to change the file path to something appropriate for your machine

 > restart
 > Export( "C:/Users/TomLeslie/Desktop/sinplot.jpg",           plot(sin(x), x=0..2*Pi, color=red)         );
 (1)
 >

## There is no concept...

of order in a mathematical set. The sets {1,2} and {2,1} are mathematically identical. The order in which the elements are displyed is meaningless.

As a "matter of convenience" Maple has to display a set in "some" order. This is often lexicographic, although I think that if the set contains certain Maple types, then the displayed order is firstly by type, then lexicographically within each type - but you shoul never, ever rely on this.

## Hard to work...

out what form of output you actually want, so the attached gives a couple of possibilities.

There are many other options you juat have to be able to specify clearly what it is that you are after

 > restart;   n := .5: N1 := 0.: N2 := 1.0: N3 := .1: lambda := .1: M := .1:   Kr := .1: Sc := 1.0: Nb := .1: Pr := 1.0: Nt := .1: R := .5:   odeSys := diff(Theta(x), x, x)+Pr*(R*diff(Theta(x), x)*f(x)+Nb*diff(Theta(x), x)*diff(Phi(x), x)+Nt*diff(Theta(x), x)^2),             N2*diff(G(x), x, x)-N1*(2*G(x)+diff(f(x), x, x))-N3*R*(diff(f(x), x)*G(x)-f(x)*diff(G(x), x)),             diff(Phi(x), x, x)+R*Sc*f(x)*diff(Phi(x), x)+Nt*diff(Theta(x), x, x)/Nb,             (1+N1)*diff(g(x), x, x)+R*(diff(g(x), x)*f(x)-g(x)*diff(f(x), x))-M*g(x)+2*Kr*diff(f(x), x),             (1+N1)*diff(f(x), x\$4)-R*(diff(f(x), x)*diff(f(x), x\$2)-f(x)*diff(f(x), x\$3))+N1*diff(G(x), x\$2)-M*diff(f(x), x\$2)-2*Kr*diff(g(x), x):   cond := f(0) = 0, (D(f))(0) = 1, g(0) = 0, Theta(0) = 1, G(0) = -n*((D@@2)(f))(0),           Phi(0) = 1, f(1) = lambda, (D(f))(1) = 0, g(1) = 0, Theta(1) = 0,           G(1) = n*((D@@2)(f))(1), Phi(1) = 0:   ans := dsolve([odeSys, cond], numeric, output=listprocedure):   eval( diff(f(x),x,x), ans)(0); # # Following two quantities are actually specifiend in cond # above as # # Theta(0) = 1,  and Phi(0) = 1 # # but you can get them from 'ans' if you want to #   eval( Theta(x), ans)(0);   eval( Phi(x), ans)(0);   plots:-odeplot( ans,                   [ [x, Theta(x)],                     [x, Phi(x)],                     [x, f(x)],                     [x, g(x)],                     [x, G(x)]                   ],                   x=0..1,                   color=[red, blue, green, black, cyan]                 );
 > interface(rtablesize=20):   Matrix( [ [x, f(x), diff(f(x), x), diff(f(x), x,x), Theta(x), Phi(x)],             seq( [ j,                    eval(f(x), ans)(j),                    eval(diff(f(x),x), ans)(j),                    eval(diff(f(x),x,x), ans)(j),                    eval(Theta(x), ans)(j),                    eval(Phi(x), ans)(j)                  ],                  j=0..1,0.1                )           ]         );   interface(rtablesize=10):
 (1)
 >

## Probably...

many ways to achieve this, of which the attached is just one

 > restart;   eqn:= g=1+c__1/x+c__2/x^2+c__3/x^3;   df_dx:= subs( op(2,eqn)=op(1,eqn),                 diff( eval( (1-1/g)*A, eqn), x)               );
 (1)
 >

## Yet another way...

- just for fun

 > restart;   with(GraphTheory):   output := table([(3, 6) = ["C-C", 0.99414575], (3, 4) = ["C-C", 1.00000000], (6, 13) = ["C-H", 1.11484141],                    (2, 7) = ["O-C", 0.99996130], (3, 5) = ["C-C", 1.60953223], (4, 7) = ["C-C", 0.99414575],                    (3, 7) = ["C-C", 1.61174619], (10, 16) = ["C-H", 0.61996904],(3, 8) = ["C-C", 0.99997800],                    (6, 7) = ["C-C", 1.60940000], (5, 12) = ["C-H", 0.61992962], (5, 6) = ["C-C", 0.99404781],                    (8, 10) = ["C-C", 0.99997800], (5, 13) = ["C-H", 0.61992962], (11, 17) = ["C-H", 0.61996904],                    (4, 5) = ["C-C", 1.60953223], (9, 15) = ["C-H", 0.62000000], (9, 11) = ["C-C", 0.99997800],                    (4, 6) = ["C-C", 1.61174619], (4, 9) = ["C-C", 0.99997800], (5, 7) = ["C-C", 0.99404781],                    (10, 11) = ["C-C", 1.00000000], (1, 6) = ["O-C", 0.99996130], (7, 12) = ["C-H", 1.11484141],                    (8, 14) = ["C-H", 0.62000000]]):   G:= Graph       ( { seq           ( [{lhs(j)[]},rhs(j)[2]],             j in remove                  ( j-> evalb(rhs(j)[1][-1]="H"),                    indices(output)=~[entries(output,'nolist')]                  )           )         }       ):   interface(rtablesize=20):   WeightMatrix(G);   interface(rtablesize=10):   DrawGraph(G, style=spring);
 >

## It is rare...

to actually see code so badly written - maybe there ought to be a prize for this!

If I remove all of the superfluous punctuation, incorrect use of commands such as sum(), product(),  unneceesary loading of the linalg() package which was deprecated about 20years ago, the unterminated "do loop" problem which Acer has already mentioned, then (to my surprise), what is left actually executes (see the attached)

Whether or not it does anything remotely useful - I have absolutely no idea

 > restart:   H := proc(z, aa, AA, bb, BB, m, n, p, q)              local mu0, i, beta0, s, Hz, v, h;              mu0 := add(BB[i], i = 1 .. q)-(add(AA[i], i = 1 .. p));              beta0 := (mul(AA[i]^(-AA[i]), i = 1 .. p))*(mul(BB[i]^(-BB[i]), i = 1 .. q));              if z <> 0 and 0 < mu0 or 0 < abs(z) and abs(z) < 1/beta0 and mu0 = 0              then s := Vector([seq(-(bb[i]-v)/BB[i], i = 1 .. m)]);                   Hz := add(add((mul(`if`(i <> h, GAMMA(bb[i]+BB[i]*s[h]), 1), i = 1 .. m))*(mul(`if`(i <> h, GAMMA(1-aa[i]-AA[i]*s[h]), 1), i = 1 .. n))*(-1)^v*z^((bb[h]+v)/BB[h])/(`if`(m < q, mul(`if`(i <> h, GAMMA(1-bb[i]-BB[i]*s[h]), 1), i = m+1 .. q), 1)*`if`(n < p, mul((`if`)(i <> h, GAMMA(aa[i]+AA[i]*s[h]), i = n+1 .. p), 1), 1)*factorial(v)*BB[h]), h = 1 .. m), v = 0 .. 10)              end if;              return Hz;         end proc:
 > H(.547, Vector([-23-j]), Vector([.5]), Vector([0, -24-j]), Vector([1, .5]), 1, 1, 1, 2);
 (1)
 >