tomleslie

13776 Reputation

19 Badges

14 years, 59 days

MaplePrimes Activity


These are answers submitted by tomleslie

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;

[b = 4, c = 6]

 

proc (z) options operator, arrow; (2*z^3+4*z^2+6)/(3*z^2+8*z) end proc

 

222.07

 

`Array(0..100, {(1) = 222.07, (2) = 147.6075359, (3) = 97.96855649, (4) = 64.87990643, (5) = 42.82682908, (6) = 28.13385322, (7) = 18.35224512, (8) = 11.85195702, (9) = 7.550115196, (10) = 4.730897051, (11) = 2.926847402, (12) = 1.840837320, (13) = 1.286750549, (14) = 1.106331044, (15) = 1.086368044, (16) = 1.086130231, (17) = 1.086130198, (18) = 1.086130198, (19) = 1.086130198, (20) = 1.086130198, (21) = 1.086130198, (22) = 1.086130198, (23) = 1.086130198, (24) = 1.086130198, (25) = 1.086130198, (26) = 1.086130198, (27) = 1.086130198, (28) = 1.086130198, (29) = 1.086130198, (30) = 1.086130198, (31) = 1.086130198, (32) = 1.086130198, (33) = 1.086130198, (34) = 1.086130198, (35) = 1.086130198, (36) = 1.086130198, (37) = 1.086130198, (38) = 1.086130198, (39) = 1.086130198, (40) = 1.086130198, (41) = 1.086130198, (42) = 1.086130198, (43) = 1.086130198, (44) = 1.086130198, (45) = 1.086130198, (46) = 1.086130198, (47) = 1.086130198, (48) = 1.086130198, (49) = 1.086130198, (50) = 1.086130198, (51) = 1.086130198, (52) = 1.086130198, (53) = 1.086130198, (54) = 1.086130198, (55) = 1.086130198, (56) = 1.086130198, (57) = 1.086130198, (58) = 1.086130198, (59) = 1.086130198, (60) = 1.086130198, (61) = 1.086130198, (62) = 1.086130198, (63) = 1.086130198, (64) = 1.086130198, (65) = 1.086130198, (66) = 1.086130198, (67) = 1.086130198, (68) = 1.086130198, (69) = 1.086130198, (70) = 1.086130198, (71) = 1.086130198, (72) = 1.086130198, (73) = 1.086130198, (74) = 1.086130198, (75) = 1.086130198, (76) = 1.086130198, (77) = 1.086130198, (78) = 1.086130198, (79) = 1.086130198, (80) = 1.086130198, (81) = 1.086130198, (82) = 1.086130198, (83) = 1.086130198, (84) = 1.086130198, (85) = 1.086130198, (86) = 1.086130198, (87) = 1.086130198, (88) = 1.086130198, (89) = 1.086130198, (90) = 1.086130198, (91) = 1.086130198, (92) = 1.086130198, (93) = 1.086130198, (94) = 1.086130198, (95) = 1.086130198, (96) = 1.086130198, (97) = 1.086130198, (98) = 1.086130198, (99) = 1.086130198, (100) = 1.086130198})`

(1)

 

Download iter.mw

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));

(alpha^2+k^2)^(1/2)*(-(alpha^2+k^2)^(1/2)*exp(k*z0)+k*exp((alpha^2+k^2)^(1/2)*z0))/((-alpha^2-k^2)*exp(k*z0)+exp((alpha^2+k^2)^(1/2)*z0)*(alpha^2+k^2)^(1/2)*k)

 

1

(1)

 

Download simplify.mw

has a size=[w,h] option.

Read the help at ?plot/options

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)
           ]
         );

  

-(1/2)*Ps*(-2*Hk*cos(theta)*sin(theta)+2*H*sin(-theta+psi)) = 0

 

 

 

Download pts.mw

 

 

 

 

In your first example,

  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!!

In your second example

 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).

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);
  

cos(x)*cosh(x)-1

 

 

[0., 4.7300407448627040260, 7.8532046240958375565, 10.995607838001670907, 14.137165491257464177, 17.278759657399481438]

(1)

 

Download rts.mw

 

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);

f(A)+f(B)

 

f(A)+f(B)+f(c)

 

f(A)+f(B)+f(c)+f(10*c)

 

sin(x)+f(A)+f(B)+(10/13)*exp(x)+cos(x)+f(c)+f(10*c)

 

f(A+B)

 

f(A+B+c)

 

f(A+B+11*c)

 

f(A+B+11*c)+sin(x)+(10/13)*exp(x)+cos(x)

(1)

 

Download addSel.mw

 

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
         );
  

diff(diff(diff(diff(diff(y(x), x), x), x), x), x) = exp(-x)*y(x)^3

 

y(0) = 1, (D(y))(0) = 1/2, ((D@@2)(y))(0) = 1/4, y(1) = exp(1/2), (D(y))(1) = (1/2)*exp(1/2)

 

 

 

Download odeProb.mw

 

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);

 

255

 

HFloat(-6.2727866497488e-17)

 

1

 

HFloat(2.0040843216548243)

 

 

 


#
# 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);

2041

 

HFloat(-0.001803176143742874)

 

1

 

HFloat(1.0)

 

 

 

 

Download autoC.mw

 

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)
        );

21356

(1)

 

Download expFile.mw

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.

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]
                );
                     

 

HFloat(-3.4606667359758356)

 

HFloat(1.0000000000000004)

 

HFloat(1.0000000000000004)

 

 

  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):

Matrix(12, 6, {(1, 1) = x, (1, 2) = f(x), (1, 3) = diff(f(x), x), (1, 4) = diff(diff(f(x), x), x), (1, 5) = Theta(x), (1, 6) = Phi(x), (2, 1) = 0, (2, 2) = 0., (2, 3) = 1.0000000000000004, (2, 4) = -3.4606667359758356, (2, 5) = 1.0000000000000004, (2, 6) = 1.0000000000000004, (3, 1) = .1, (3, 2) = 0.835702022307008e-1, (3, 3) = .6800631923030176, (3, 4) = -2.9410234948062945, (3, 5) = .9057823964917024, (3, 6) = .8909960606660462, (4, 1) = .2, (4, 2) = .13771757959262818, (4, 3) = .4112908898410515, (4, 4) = -2.4365617798889927, (4, 5) = .8100219534215602, (4, 6) = .7839719522758213, (5, 1) = .3, (5, 2) = .1674902374939747, (5, 3) = .19238609178184773, (5, 4) = -1.9430829821978786, (5, 5) = .7129641055005085, (5, 6) = .6789638114289943, (6, 1) = .4, (6, 2) = .17782543711410675, (6, 3) = 0.22408618058967926e-1, (6, 4) = -1.4576065601900232, (6, 5) = .6147473554036548, (6, 6) = .5759767479614302, (7, 1) = .5, (7, 2) = .17357961237278668, (7, 3) = -0.9933396660786277e-1, (7, 4) = -.9781164239030892, (7, 5) = .5154199662489422, (7, 6) = .4750008212147106, (8, 1) = .6, (8, 2) = .15954873659849525, (8, 3) = -.1733704061714211, (8, 4) = -.503314426530307, (8, 5) = .41495891428165227, (8, 6) = .37602233302662325, (9, 1) = .7, (9, 2) = .14048150220512115, (9, 3) = -.20012556462713801, (9, 4) = -0.323820966655486e-1, (9, 5) = .31329030919518824, (9, 6) = .2790312873754467, (10, 1) = .8, (10, 2) = .12108769733119493, (10, 3) = -.17995681820337722, (10, 4) = .4352496015845534, (10, 5) = .21031099360573124, (10, 6) = .18402544968397186, (11, 1) = .9, (11, 2) = .10604408233364852, (11, 3) = -.11316778271460856, (11, 4) = .9001235496331798, (11, 5) = .10591140472049596, (11, 6) = 0.9101106171886827e-1, (12, 1) = 1.0, (12, 2) = 0.9999999999999996e-1, (12, 3) = 0., (12, 4) = 1.362974267178911, (12, 5) = 0., (12, 6) = 0.})

(1)

 

Download odeSol.mw

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)
              );

g = 1+c__1/x+c__2/x^2+c__3/x^3

 

(-c__1/x^2-2*c__2/x^3-3*c__3/x^4)*A/g^2

(1)

 

Download oddSub.mw

 - 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);

Matrix(11, 11, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = .99996130, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (1, 10) = 0, (1, 11) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = .99996130, (2, 8) = 0, (2, 9) = 0, (2, 10) = 0, (2, 11) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 1.00000000, (3, 5) = 1.60953223, (3, 6) = .99414575, (3, 7) = 1.61174619, (3, 8) = .99997800, (3, 9) = 0, (3, 10) = 0, (3, 11) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 1.00000000, (4, 4) = 0, (4, 5) = 1.60953223, (4, 6) = 1.61174619, (4, 7) = .99414575, (4, 8) = 0, (4, 9) = .99997800, (4, 10) = 0, (4, 11) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 1.60953223, (5, 4) = 1.60953223, (5, 5) = 0, (5, 6) = .99404781, (5, 7) = .99404781, (5, 8) = 0, (5, 9) = 0, (5, 10) = 0, (5, 11) = 0, (6, 1) = .99996130, (6, 2) = 0, (6, 3) = .99414575, (6, 4) = 1.61174619, (6, 5) = .99404781, (6, 6) = 0, (6, 7) = 1.60940000, (6, 8) = 0, (6, 9) = 0, (6, 10) = 0, (6, 11) = 0, (7, 1) = 0, (7, 2) = .99996130, (7, 3) = 1.61174619, (7, 4) = .99414575, (7, 5) = .99404781, (7, 6) = 1.60940000, (7, 7) = 0, (7, 8) = 0, (7, 9) = 0, (7, 10) = 0, (7, 11) = 0, (8, 1) = 0, (8, 2) = 0, (8, 3) = .99997800, (8, 4) = 0, (8, 5) = 0, (8, 6) = 0, (8, 7) = 0, (8, 8) = 0, (8, 9) = 0, (8, 10) = .99997800, (8, 11) = 0, (9, 1) = 0, (9, 2) = 0, (9, 3) = 0, (9, 4) = .99997800, (9, 5) = 0, (9, 6) = 0, (9, 7) = 0, (9, 8) = 0, (9, 9) = 0, (9, 10) = 0, (9, 11) = .99997800, (10, 1) = 0, (10, 2) = 0, (10, 3) = 0, (10, 4) = 0, (10, 5) = 0, (10, 6) = 0, (10, 7) = 0, (10, 8) = .99997800, (10, 9) = 0, (10, 10) = 0, (10, 11) = 1.00000000, (11, 1) = 0, (11, 2) = 0, (11, 3) = 0, (11, 4) = 0, (11, 5) = 0, (11, 6) = 0, (11, 7) = 0, (11, 8) = 0, (11, 9) = .99997800, (11, 10) = 1.00000000, (11, 11) = 0})

 

 

 

Download grph.mw

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);

.578683260/GAMMA(-.5*v+25+j)

(1)

 

Download reallyBadCode.mw

 

5 6 7 8 9 10 11 Last Page 7 of 207