Advanced Search

tomleslie

13821 Reputation

20 Badges

14 years, 244 days
Contact tomleslie

MaplePrimes Activity


These are answers submitted by tomleslie

Something like...

tomleslie 13821
February 21 2023
0 0

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

Well...

tomleslie 13821
February 21 2023
2 0

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

The plot() command...

tomleslie 13821
February 19 2023
0 0

has a size=[w,h] option.

Read the help at ?plot/options

After fixing...

tomleslie 13821
February 19 2023
0 0

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

 

 

 

 

Some observations...

tomleslie 13821
February 17 2023
0 1

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

More roots...

tomleslie 13821
February 16 2023
0 0

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

 

It is not too difficult...

tomleslie 13821
February 15 2023
0 0

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

 

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

tomleslie 13821
February 14 2023
0 0

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

 

Interesting one...

tomleslie 13821
February 14 2023
2 4

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

 

Trivial...

tomleslie 13821
February 12 2023
2 0

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

There is no concept...

tomleslie 13821
February 10 2023
0 3

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

tomleslie 13821
February 10 2023
0 1

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

Probably...

tomleslie 13821
February 09 2023
2 0

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

Yet another way...

tomleslie 13821
February 09 2023
1 3

 - 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

It is rare...

tomleslie 13821
February 08 2023
1 3

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