tomleslie

6237 Reputation

17 Badges

10 years, 105 days

MaplePrimes Activity


These are answers submitted by tomleslie

cos (s)he has posted a correct answer.

Unfortunatley for you it is unliklet that anyoen is going to download a package (Gym-pakke ) to do simple linear fits, becuase

  1. it's in Danish - for me English is fine, French is OK and German, maybe - but Danish? Sorry it isn't going to happen
  2. who needs some add-on packag to do a simple linear fit? This is trivial to do with standard Maple commands as mmcdara has already shown

See also the attached

  restart:
  with(plots):
#
# The data
#
  X:=<0,    10,  20,  30,  40,  50,  60,  70>:
  Y:=<690, 716, 742, 770, 798, 825, 850, 877>:
#
# The linear fit
#
  f:=Statistics:-LinearFit(a+b*x, X, Y, x);
#
# Plots of the data and the 'fit'
#
  plots:-display( [ pointplot(X, Y, symbol=solidcircle, symbolsize=20, color=red),
                    plot(f, x=X[1]..X[-1])
                  ]
                );

HFloat(689.4999999999998)+HFloat(2.6857142857142895)*x

 

 

 

 

 

Download linFit.mw

was the use of the parameter 'a', for which there is no numerical value - so I highlighted its use in the attached and (arbitrarliy) gave it the value 1.

I fixed a few other typos

restart

with(plots)

inf := 50

lambda := L(t)+k*A(t)

odes := diff(S(t), t) = Lambda-(`&beta;__1`+`&beta;__2`)*S(t)*lambda-mu*S(t)+r*R(t), diff(E(t), t) = `&beta;__1`*S(t)*lambda-(alpha+mu+d__1+(1-p)*epsilon+p*eta)*E(t), diff(L(t), t) = p*eta*E(t)+e*A(t)-(`&sigma;__2`+mu+d__3+y__2)*L(t), diff(A(t), t) = (1-p)*epsilon*E(t)-(`&sigma;__1`+mu+d__2+y__1+e)*A(t), diff(R(t), t) = y__2*L(t)+y__1*A(t)-(mu+r)*R(t), diff(V(t), t) = `&beta;__2`*S(t)*lambda+alpha*E(t)+`&sigma;__1`*A(t)+`&sigma;__2`*L(t)-a*V(t)

ICs := S(0) = S__0, E(0) = E__0, L(0) = L__0, A(0) = A__0, R(0) = R__0, V(0) = V__0

parameters := [S__0 = 1000000, E__0 = 1, L__0 = 0, A__0 = 0, R__0 = 0, V__0 = 0, Lambda = .22, mu = .22, `&beta;__1` = 0.3e-1, `&beta;__2` = 0.3e-1, r = .5, alpha = .6, p = .1, k = 0.4e-1, y__1 = .8, y__2 = .8, d__1 = 0.2e-1, d__2 = 0.6e-1, d__3 = .2, `&sigma;__1` = 0.2e-1, `&sigma;__2` = 0.2e-1, epsilon = 0.2e-1, eta = 0.2e-1, e = 0.1e-1, a = 1]

S1 := dsolve(eval([odes, ICs], parameters), numeric)

P1 := odeplot(S1, [[t, S(t)]], 0 .. inf, linestyle = 1, color = blue, labels = ["time", "Population"], axes = "boxed"); P2 := odeplot(S1, [[t, E(t)]], 0 .. inf, linestyle = 2, color = brown, labels = ["time", "Population"], axes = "boxed"); P3 := odeplot(S1, [[t, L(t)]], 0 .. inf, linestyle = 3, color = red, labels = ["time", "Population"], axes = "boxed"); P4 := odeplot(S1, [[t, A(t)]], 0 .. inf, linestyle = 4, color = green, labels = ["time", "Population"], axes = "boxed"); P5 := odeplot(S1, [[t, R(t)]], 0 .. inf, linestyle = 5, color = purple, labels = ["time", "Population"], axes = "boxed")

 

 

 

 

 

``


 

Download anotherCOVID.mw

Or at least something very like it?

The attached performs the calculations and animations - alhtough combining the pre-launch and post-launch animations is a bit uninformative becuase of the large discrepancy in both space and time in the two regions

  restart;
  with(plots):
  params:= [alpha=Pi/6, l=40, V__B=4.5, d=3, g=9.81]:
#
# Pre-launch solution and animation
#
  ode1:= diff(x__1(t),t,t)=P-m*g*sin(alpha):
  ics:=  x__1(0)=0, D(x__1)(0)=0:
  sol1:=dsolve( [ ode1, ics ] );
  s1:=eval( sol1, [ t=tau, x__1(t)=l]);
  s2:=eval( diff( sol1, t ), [ diff(x__1(t),t)=V__B, t=tau])/tau:
  sol2:=algsubs( rhs(s2)=lhs(s2), [s1, sol1]);
  tau__val:=eval(isolate( sol2[1], tau), params);
  sol3:=eval(sol2[2], [ params[], tau__val]):
  animate( pointplot,
           [ eval([(rhs(sol3)-l)*sin(alpha),(rhs(sol3)-l)*cos(alpha)], params),
             symbol=solidcircle,
             symbolsize=20
           ],
           t=0..rhs(tau__val),
           frames=20,
           trace=20
        );

x__1(t) = (1/2)*(P-m*g*sin(alpha))*t^2

 

l = (1/2)*(P-m*g*sin(alpha))*tau^2

 

[l = (1/2)*tau*V__B, x__1(t) = (1/2)*t^2*V__B/tau]

 

tau = 17.77777778

 

 

#
# Post-launch solutioon and animation
#
  ode2:= diff(x(t),t,t)=0,
         diff(y(t),t,t)=-g:
  ics:= x(0)=0,
        y(0)=0,
        D(x)(0)=V__B*cos(alpha),
        D(y)(0)=V__B*sin(alpha):
  sol4:= dsolve( [ ode2, ics ] ):
  sol5:=eval(eliminate( eval( sol4, [ x(t)=d, y(t)=-h, t=T] ), T), params):
  h__val=evalf(rhs(isolate( sol5[2][], h)));
  animate( pointplot,
           [ eval([rhs~(sol4)[]], params),
             symbol=solidcircle,
             symbolsize=20
           ],
           t=0..rhs(sol5[1][]),
           frames=20,
           trace=20
        );

h__val = 1.174615859

 

 

#
# It is possible to combine the two animations, but
#
# 1. the vertical scales are so different, (~40 pre-
#    launch, and 1 post-launch), and
# 2. the duration of the pre-launch and post-launch
#    periods are so different ( ~18secs pre-launch
#    and ~0.8secs post launch
#
# means that the combined animation doesn't "look"
# very good
#
  animate
  ( pointplot,
    [ piecewise
      ( t<rhs(tau__val),
        eval([(rhs(sol3)-l)*sin(alpha),(rhs(sol3)-l)*cos(alpha)], params),
        eval(subs( t=t-rhs(tau__val), eval([rhs~(sol4)[]], params)), params)
      ),
      symbol=solidcircle,
      symbolsize=20
    ],
    t=0..rhs(tau__val)+rhs(sol5[1][]),
    frames=100,
    trace=100
  );

 

 

Download launch.mw

 

in the Mapleprimes toolbar to upload your worksheets

As you presented it I had to do a lot of "hacking" just to figure out what tou *might* want - and I *may* have got it wrong.

So all I can say is that the attached worksheet *might* perform the calculation(s) you want!!

  restart;
  with(plots):
  hz2 := phi1(t)+arcsin((h+r1*sin(phi1(t)))/r2);
  deq := diff(phi1(t), t, t)-hz2 = 0;
#
# Two sets of initial conditions
#
  ics1 := phi1(0) = -arcsin(h/r1), (D(phi1))(0) = 0;
  ics2 := phi1(0) = 0, D(phi1)(0) = -0.63;

phi1(t)+arcsin((h+r1*sin(phi1(t)))/r2)

 

diff(diff(phi1(t), t), t)-phi1(t)-arcsin((h+r1*sin(phi1(t)))/r2) = 0

 

phi1(0) = -arcsin(h/r1), (D(phi1))(0) = 0

 

phi1(0) = 0, (D(phi1))(0) = -.63

(1)

#
# Solve/plot ode with initial conditions ic1
#
  sol1 := dsolve
          ( eval
            ( [ deq, ics1 ],
              [ r1 = 8, r2 = 8, r3 = 1, h = 5,
                m2 = 1, m3 = 20, theta3 = 20
              ]
            ),
            phi1(t),
            numeric
          ):
   odeplot( sol1,
            [ t, phi1(t) ],
            t=0..2
          );

 

#
# Solve/plot ode with initial conditions ic2
#
   sol2:= dsolve
          ( eval
            ( [ deq, ics2 ],
              [ r1 = 8, r2 = 8, r3 = 1, h = 5,
                m2 = 1, m3 = 20, theta3 = 20
              ]
            ),
            phi1(t),
            numeric
          ):
   odeplot( sol2,
            [ t, phi1(t) ],
            t=0..2
          );

 

 

Download solveODE.mw

Using the DETools:-DEplot() may be a bit "overkill" if you just want the solution curves. See the attached for an alternative

  restart;
  ode1 := diff(P(t), t) = 0.2*P(t) - 300:
  sol1:=dsolve([ ode1, P(0)=1300]);
  sol2:=dsolve([ ode1, P(0)=1800]);
  plot(rhs~([sol1, sol2]), t=0..20);

P(t) = 1500-200*exp((1/5)*t)

 

P(t) = 1500+300*exp((1/5)*t)

 

 

s

s

(1)

 


 

Download odeplt.mw

to figure out precisely what a poster needs when they describe a plot - but maybe as shown in the attached?

BTW - the "gridlines" will not appear in an actual Maple worksheet (unless you set the option gridlines=true). I have given up wonderiing why they always appear in any plots in worksheets uploaded to this site

  l1:= plottools:-line( [-Pi/4,0], [-Pi/4, sec(-Pi/4)], color=black, thickness=5):
  l2:= plottools:-line( [ Pi/4,0], [ Pi/4, sec( Pi/4)], color=black, thickness=5):
  p:=  plots:-shadebetween(sec(x),0,x=-Pi/4..Pi/4, color=red):
  plots:-display([l1,l2,p], view=[-Pi/3..Pi/3, 0..1.1*sec(Pi/4)]);

 

 


 

Download secplt.mw

The attached is a seriously "butchered" version of your code. I removed stuff wihch wasn't doing anything, combined some loops for efficiency, etc, etc - but the resulting version is still embarrassingly bad.

However it will perform the required calculations for any supplied value of the parameter 'alpha', returning a "2D" graph (technically a 3D spacecurve) and a 3D graph.

These spacecurves and surface plots can then be combined in any way you want. Several possibilities are shown in the attached.

  restart:
  getPlots:= proc( alpha1, col1, col2)
                   local k:=2, M:=2, K:=2^(k-1),
                         N:=K*M, beta:=1, i, X, T, r, Phi_mxm,
                         key:=4, rho:=1, mmm:=1, sigma:=1,
                         lambda:=(m,alpha1,beta)->sqrt((2*m+alpha1+beta+1)*GAMMA(2*m+alpha1+beta+1)*m!/(2^(alpha1+beta+1)*GAMMA(m+alpha1+1)*GAMMA(m+beta+1))),
                         psii:=(n,m,x)-> piecewise( (n-1)/K <= x and x <= n/K,
                                                     2^(k/2)*lambda(m,alpha1,beta)*simplify(JacobiP(m,alpha1,beta,2^(k)*x-2*n+1)),
                                                     0
                                                  ),
                         w:=(n,x)->(1+x)^(1/gama-1),
                         omega:=(n,x)->w(n,2^(k)*x-2*n+1),
                         psi:=t->Matrix(N,1,[seq(seq(psii(i,j,t),j=0..M-1),i=1...K)] ),
                         PP:=alpha-> Phi_mxm.P(k,M,alpha).MatrixInverse(Phi_mxm),
                         u_exact:=(x,t)-> 1/(1+exp(sqrt(key/6)*x-(5/6)*key*t))^2,
                         f1:=unapply(u_exact(x,0),x),
                         g1:=unapply(u_exact(0,t),t),
                         g2:=unapply(u_exact(1,t),t),
                         U:=Matrix( N,N,symbol=u),
                         wwxx, wwt,ww,MC, PDEson, j, coz, UU,
                         P:= proc( k,M,nn )
                                   local PB,m,i,p,j,xi;
                                   m:=M*(2^(k-1)):
                                   xi:=(i,n)->((i+1)^(n+1)-2*i^(n+1)+(i-1)^(n+1));
                                   PB:=Matrix(m,m):
                                   for i from 1 to m do
                                       p:=0:
                                       for j from 1 to m do
                                           if   i=j
                                           then PB[i,j]:=1;
                                           fi:
                                           if   i<j
                                           then p:=p+1:
                                                PB[i,j]:= xi(p,nn);
                                           fi:
                                       end do;
                                       p:=1:
                                   end do:
                                   PB:=1/m^nn/GAMMA(nn+2)*PB;
                                   return PB;
                             end proc:
                   uses plots, LinearAlgebra;
                   for i from 1 to N do
                       X[i]:=evalf((2*i-1)/((2^k)*M)):
                       T[i]:=evalf((2*i-1)/((2^k)*M)):
                       r[i]:=evalf(psi(T[i])):
                   end do:

                   Phi_mxm:=Matrix([seq(r[i],i=1...N)]);
                   wwxx:= diff(f1(x),x,x)+(psi(x)^+.U.PP(alpha1).psi(t))(1):
                   wwt:= diff(g1(t),t)+x*(diff(g2(t),t)-diff(g1(t),t)-(psi(1)^+.PP(2)^+.U.psi(t))+psi(x)^+.PP(2)^+.U.psi(t))(1) :
                   ww:= f1(x)+g1(t)-g1(0)+x*(g2(t)-g2(0)-g1(t)+g1(0)-(psi(1)^+.PP(2)^+.U.PP(alpha1).psi(t)))(1)+(psi(x)^+.PP(2)^+.U.PP(alpha1).psi(t))(1) :
                   MC:=unapply(wwt -rho*wwxx -key*ww.(1-(1/mmm)*ww^sigma ),x,t):
                   PDEson:=Matrix(N, N, (i,j)-> simplify( evalf( MC(X[i],T[j] ) ))=0):
                   coz:=fsolve( { seq( seq(PDEson(i,j),i=1..N ),j=1..N )}):
                   U:=subs(op(coz),U):
                   UU:=(x,t)->evalf(f1(x)+g1(t)-g1(0)+x*(g2(t)-g2(0)-g1(t)+g1(0)-(psi(1)^+.PP(2)^+.U.PP(alpha1).psi(t)))(1)+(psi(x)^+.PP(2)^+.U.PP(alpha1).psi(t))(1)):
                   return plot3d(UU(x,t),x=0..1,t=0..1,color=col1, transparency=0.5),
                          plots:-spacecurve([0.5, t, UU(0.5,t)],t=0..1, color=col2, thickness=3):
             end proc:

  p1,p2:=getPlots(1, red, blue):
  p3,p4:=getPlots(0.3, green, orange):
  plots:-display([p1,p2]);
  plots:-display([p3,p4]);
  plots:-display([p1,p3]);
  plots:-display([p2,p4]);

 

 

 

 

 


 

Download grph23v2.mw

 

as shown in the final group of the attached


 

 

restart:
with(orthopoly):
with(LinearAlgebra):
with(plots):
interface(rtablesize=20):

k:=2:
M:=2:
K:=2^(k-1):
N:=K*M:

 

 

alpha:=1:


alpha1:=alpha: 
beta:=1:
lambda:=(m,alpha1,beta)->sqrt((2*m+alpha1+beta+1)*GAMMA(2*m+alpha1+beta+1)*m!/(2^(alpha1+beta+1)*GAMMA(m+alpha1+1)*GAMMA(m+beta+1)));
psii:=(n,m,x)->piecewise((n-1)/K <= x and x <= n/K,  2^(k/2) *lambda(m,alpha1,beta)*simplify(JacobiP(m,alpha1,beta,2^(k)*x-2*n+1)), 0);
local i,j: psi:=(x)->Array([seq(seq(psii(i,j,x),j=0..M-1),i=1..K)] ):

w:=(n,x)->(1+x)^(1/gama-1):
omega:=(n,x)->w(n,2^(k)*x-2*n+1):

lambda := proc (m, alpha1, beta) options operator, arrow; sqrt((2*m+alpha1+beta+1)*GAMMA(2*m+alpha1+beta+1)*factorial(m)/(2^(alpha1+beta+1)*GAMMA(m+alpha1+1)*GAMMA(m+beta+1))) end proc

 

proc (n, m, x) options operator, arrow; piecewise((n-1)/K <= x and x <= n/K, 2^((1/2)*k)*lambda(m, alpha1, beta)*simplify(JacobiP(m, alpha1, beta, 2^k*x-2*n+1)), 0) end proc

(1.1)


psi:=t->Matrix(N,1,[seq(seq(psii(i,j,t),j=0..M-1),i=1...K)] ):
 
 


for i from 1 to N do
X[i]:=evalf((2*i-1)/((2^k)*M)):
end do:

for i from 1 to N do
T[i]:=evalf((2*i-1)/((2^k)*M)):
end do:

for i from 1 to N do
r[i]:=evalf(psi(T[i])):
end do:
Phi_mxm:=Matrix([seq(r[i],i=1...N)]);
 

Matrix(4, 4, {(1, 1) = 1.732050808, (1, 2) = 1.732050808, (1, 3) = 0., (1, 4) = 0., (2, 1) = -3.872983346, (2, 2) = 3.872983346, (2, 3) = 0., (2, 4) = 0., (3, 1) = 0., (3, 2) = 0., (3, 3) = 1.732050808, (3, 4) = 1.732050808, (4, 1) = 0., (4, 2) = 0., (4, 3) = -3.872983346, (4, 4) = 3.872983346})

(1.2)

 

P:=proc(k,M,nn) local PB,m,i,p,j,xi;
m:=M*(2^(k-1)):
xi:=(i,n)->((i+1)^(n+1)-2*i^(n+1)+(i-1)^(n+1));
PB:=Matrix(m,m):
for i from 1 to m do
p:=0:
for j from 1 to m do
if i=j then PB[i,j]:=1;
 fi:
if i<j then p:=p+1:
PB[i,j]:= xi(p,nn);
 fi:
end do;
p:=1:
end do:
PB:=1/m^nn/GAMMA(nn+2)*PB;

return PB;
end proc:
PP:=alpha->Phi_mxm.P(k,M,alpha).MatrixInverse(Phi_mxm);

proc (alpha) options operator, arrow; `.`(Phi_mxm, P(k, M, alpha), LinearAlgebra:-MatrixInverse(Phi_mxm)) end proc

(1.3)


key:=4:
rho:=1:
mmm:=1:
sigma:=1:

u_exact:=(x,t)-> 1/(1+exp(sqrt(key/6)*x-(5/6)*key*t))^2 ;
f1:=unapply(u_exact(x,0),x):
g1:=unapply(u_exact(0,t),t);
g2:=unapply(u_exact(1,t),t);

 

U:=Matrix( N,N,symbol=u);
wwxx:= diff(f1(x),x,x)+(psi(x)^+.U.PP(alpha).psi(t))(1):
wwt:= diff(g1(t),t)+x*(diff(g2(t),t)-diff(g1(t),t)-(psi(1)^+.PP(2)^+.U.psi(t))+psi(x)^+.PP(2)^+.U.psi(t))(1) :
ww:= f1(x)+g1(t)-g1(0)+x*(g2(t)-g2(0)-g1(t)+g1(0)-(psi(1)^+.PP(2)^+.U.PP(alpha).psi(t)))(1)+(psi(x)^+.PP(2)^+.U.PP(alpha).psi(t))(1) :

MC:=unapply(wwt -rho*wwxx -key*ww .(1-(1/mmm)*ww^sigma ),x,t):

u_exact := proc (x, t) options operator, arrow; 1/(1+exp(sqrt((1/6)*key)*x-(5/6)*key*t))^2 end proc

 

g1 := proc (t) options operator, arrow; 1/(1+exp(-(10/3)*t))^2 end proc

 

g2 := proc (t) options operator, arrow; 1/(1+exp((1/3)*sqrt(6)-(10/3)*t))^2 end proc

 

Matrix(%id = 18446744074368602590)

(1)

PDEson:=Matrix(N,N):

for i from 1 to N do
for j from 1 to N do
PDEson(i,j):=simplify(
evalf(

MC(X[i],T[j])


)
)=0:
end do:
end do:

coz:=fsolve({seq(seq(PDEson(i,j),i=1..N ),j=1..N )}):
U:=subs(op(coz),U):
 


UU:=(x,t)->evalf(f1(x)+g1(t)-g1(0)+x*(g2(t)-g2(0)-g1(t)+g1(0)-(psi(1)^+.PP(2)^+.U.PP(alpha).psi(t)))(1)+(psi(x)^+.PP(2)^+.U.PP(alpha).psi(t))(1)):


with(plots):
graphic_3D:=plot3d(UU(x,t),x=0..1,t=0..1,color=gray, transparency=0.5);

 

graphic_2D:=plots:-spacecurve([0.5, t, UU(0.5,t)],t=0..1, color=green, thickness=3):
plots:-display([graphic_2D, graphic_3D]);

 

 

 

 

 


 

Download grph23.mw

  1. either use evalc() which wii essentially simplify the expression based on the assumption that all unknown names are 'real, or
  2. use Maple's assume() command, to specify that both 'x' and 'lambda' are real

Both mehod's are shown in the attached

restart;
conjugate(1-I*x*lambda^2)+1-I*x*lambda^2;
evalc(%);

2+I*conjugate(x*lambda^2)-I*x*lambda^2

 

2

(1)

restart:
interface(showassumed=0):
assume(lambda::real, x::real);
conjugate(1-I*x*lambda^2)+1-I*x*lambda^2;

2

(2)

 

Download comp.mw

 

is worth reading!!!

Excerpts from the help at ?fsolve/details states

  • For a single polynomial equation of one variable, the fsolve command computes all real (non-complex) roots.

  • For a general equation or system of equations or procedures, the fsolve command computes a single real root.

Which part of these statements don't you understand?

 

 

You can apply a function f() to a list 'L' of values, by using the elementwise operator '~', as in  f~(L).

See the attached for a couple of ways

  restart:
  fprime_expr:=x^sin(x)*(cos(x)*ln(x)+sin(x)/x);
  RootFinding:-Analytic(diff(fprime_expr,x,x)=0, re=0..10, im=-1..1);
  L:=sort(select(type,[%],realcons));

  f:=x->x^sin(x);
#
# Apply 'f' elementwise to L
#
  f~(L);

x^sin(x)*(cos(x)*ln(x)+sin(x)/x)

 

4.78206109907178, 3.48122397295652, 2.16436178885414, .748344075670115-.240667154995574*I, .748344075670115+.240667154995574*I, 9.04887132811185, 6.77657190535075, 7.93089684792655

 

[2.16436178885414, 3.48122397295652, 4.78206109907178, 6.77657190535075, 7.93089684792655, 9.04887132811185]

 

proc (x) options operator, arrow; x^sin(x) end proc

 

[1.896584783, .6599753075, .2099102798, 2.475003072, 7.882490202, 2.244817923]

(1)

  restart:
  fprime_expr:=x^sin(x)*(cos(x)*ln(x)+sin(x)/x):
  getRoots:=proc( g::procedure)
                  local x__old:=0,
                        x__new:=[],
                        rt:
                  while true do
                        rt:= RootFinding:-NextZero(g, x__old):
                        if   rt <10
                        then x__new:=[x__new[], rt];
                             x__old:=x__new[-1];
                        else break;
                        fi;
                  od:
                  return x__new
            end proc:
  L:=getRoots( unapply( diff(fprime_expr,x,x), x) );
  f:=x->x^sin(x);
#
# Apply 'f' elementwise to L
#
  f~(L);   

[2.164361788, 3.481223972, 4.782061099, 6.776571905, 7.930896847, 9.048871328]

 

proc (x) options operator, arrow; x^sin(x) end proc

 

[1.896584783, .6599753083, .2099102798, 2.475003072, 7.882490202, 2.244817923]

(2)

 


 

Download fList.mw

  1. Either use' evalc()', so that Maple is forced to evaluate expressions with the assumption that al unknown names are real
  2. Or, use the 'assume()' facility to explicitly declare certain variables as 'real'

Both methods are shown in the attached

#
# Use 'evalc()' to ensure that all expressions
# evaluate with the 'assumption' that any anknown
# names are real
#
  restart:
  alias(conj=conjugate);

  d := a*x+b*y-c = 0;
  z := x+I*y;
  evalc(z+conj(z));
  evalc(z-conj(z));

  d := evalc((1/2)*a*(z+conj(z))+b*(z-conj(z))/(2*I)-c);

  is(d = (1/2)*evalc(z*(a-I*b)+conj(z)*(a+I*b)-2*c));

  varpi:= a+I*b;
  is(d = (1/2)*evalc(z*conj(varpi)+conj(z)*varpi-2*c));

conj

 

a*x+b*y-c = 0

 

x+I*y

 

2*x

 

(2*I)*y

 

a*x+b*y-c

 

true

 

a+I*b

 

true

(1)

#
# Use assumptions to specify certain names are
# real
#
  restart:
  alias(conj=conjugate):
  assume(a::real, b::real, x::real, y::real):
  interface(showassumed=0);

  d:= a*x+b*y-c = 0;
  z:= x+I*y;
  z+conj(z);
  z-conj(z);

  d:= (1/2)*a*(z+conj(z))+b*(z-conj(z))/(2*I)-c;

  is(d = (1/2)*(z*(a-I*b)+conj(z)*(a+I*b)-2*c));

  varpi:= a+I*b;
  is(d = (1/2)*(z*conj(varpi)+conj(z)*varpi-2*c));

1

 

a*x+b*y-c = 0

 

x+I*y

 

2*x

 

(2*I)*y

 

a*x+b*y-c

 

true

 

a+I*b

 

true

(2)

 

Download cplex.mw

For a finite number of terms it is generally a better idea to use add() rather than sum() - particularly if there is absolutely no prospect that the call to sum() wii return a "closed form" expression.

In principle, if sum() cannot find a "closed-form" solution, the it should default to add() - but there seems to be a difference: see the attached.

So I don;t think this is actually a 0^0 issue

  restart;
  f:=x->exp(x);
  p:=n->(x->sum( a[j]*x^j, j=0..n));
  p2:=p(2);
  bed:={p2(0)=f(0), D(p2)(0)=D(f)(0), (D@@2)(p2)(0)=(D@@2)(f)(0)}

proc (x) options operator, arrow; exp(x) end proc

 

proc (n) options operator, arrow; proc (x) options operator, arrow; sum(a[j]*x^j, j = 0 .. n) end proc end proc

 

proc (x) options operator, arrow; sum(a[j]*x^j, j = 0 .. 2) end proc

 

{0 = 1, a[1] = 1, 2*a[2] = 1}

(1)

  restart;
  f:=x->exp(x);
  p:=n->(x->add( a[j]*x^j, j=0..n));
  p2:=p(2);
  bed:={p2(0)=f(0), D(p2)(0)=D(f)(0), (D@@2)(p2)(0)=(D@@2)(f)(0)}

proc (x) options operator, arrow; exp(x) end proc

 

proc (n) options operator, arrow; proc (x) local j; options operator, arrow; add(a[j]*x^j, j = 0 .. n) end proc end proc

 

proc (x) local j; options operator, arrow; add(a[j]*x^j, j = 0 .. 2) end proc

 

{a[0] = 1, a[1] = 1, 2*a[2] = 1}

(2)

  

 

Download sumadd.mw

"inflection points are where the gradient of the function is zero - that's the first derivative, not the second as given by Kitonum.

The attached worksheet plots

  1. the function fprime_expr
  2. the "roots" of the function fprime_expr, ie the points at which fprime_expr=0
  3. the "inflection points" of the function fprime_expr, ie the points at which diff(fprime_expr, x)=0

  restart:
  with(plots):
  fprime_expr:=x^sin(x)*(cos(x)*ln(x)+sin(x)/x):
  getRoots:=proc( g::procedure)
                  local x__old:=0,
                        x__new:=[],
                        rt:
                  while true do
                        rt:= RootFinding:-NextZero(g, x__old):
                        if   rt <10
                        then x__new:=[x__new[], rt];
                             x__old:=x__new[-1];
                        else break;
                        fi;
                  od:
                  return x__new
            end proc:       

#
# Plot the curve given by fprime_expr
#
  p1:= plot(fprime_expr, x=0..10):
#
# Plot the "roots" of the function fprime_expr
# ie the points where fprime_expr=0
#
  p2:= pointplot( [ seq( [j, eval(fprime_expr, x=j)],
                        j in getRoots( unapply(fprime_expr,x))
                      )
                  ],
                  symbol=solidcircle,
                  symbolsize=20,
                  color=red
                ):

#
# Plot the "inflection points" of the function fprime_expr
# ie the points where diff( fprime_expr, x) = 0
#
  p3:= pointplot( [ seq( [j, eval(fprime_expr, x=j)],
                         j in getRoots( unapply(diff(fprime_expr,x),x))
                       )
                  ],
                  symbol=solidcircle,
                  symbolsize=20,
                  color=blue
                ):
  display([p1,p2, p3]);

 

 


 

Download rootInf.mw

 

 to produce a simple animation!

I mean what is wrong with the attached??

  plots:-animate(plot, [[3*t, 4*t^2+1, t=0..tau]], tau=0..Pi);

 

 


 

Download animplt.mw

 

1 2 3 4 5 6 7 Last Page 1 of 131