Advanced Search

dharr

Dr. David Harrington

8235 Reputation

22 Badges

20 years, 340 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
Contact Dr. David Harrington
I am a retired professor of chemistry at the University of Victoria, BC, Canada. My research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

MaplePrimes Activity


These are answers submitted by dharr

poles in complex plane?...

dharr 8235
November 17 2020
0 0

You are assuming we know what field you are working in, but I don't really know what you want. A simple plot shows that the range of g is -infinity..infinity, so that is the answer to the direct question asked. But now you mention evalc, so perhaps you want the poles in the complex plane?
 

> 

g:=1/(x+x^2+x^3+2)

1/(x^3+x^2+x+2)

> 

plot(g,x=-10..10);

Poles of g in complex plane

> 

fsolve(1/g,complex);
plots:-complexplot([%],style=point,symbolsize=15,symbol=solidcircle,color=red);

HFloat(-1.3532099641993245), HFloat(0.17660498209966224)-HFloat(1.202820819285479)*I, HFloat(0.17660498209966224)+HFloat(1.202820819285479)*I

> 

 

 

Download gain.mw

pencil of conics through four points...

dharr 8235
November 13 2020
1 2

Here is how I would do it - can be refined for nicer spacing etc.

> 

restart;
> 

with(plots):local D;

Four points. We want a family of conics passing through these four points

> 

pts:=[[1,3/2],[-1/2,1],[4/3,-3/4],[-1,-1]];
p1:=plot(pts,style=point,colour=red,symbol=solidcircle,symbolsize=15):display(p1);

[[1, 3/2], [-1/2, 1], [4/3, -3/4], [-1, -1]]

General Conic, scaled so constant (F) = 1

> 

gen:=A*x^2+B*x*y+C*y^2+D*x+E*y+1=0;

A*x^2+B*x*y+C*y^2+D*x+E*y+1 = 0

Equations to be satisfied at the four points and at a general 5th point (X,Y), which we will vary

> 

eqns:=[eval(gen,[x=X,y=Y]),map(i->eval(gen,{x=i[1],y=i[2]}),pts)[]];
solve(eqns,[A,B,C,D,E]);
con:=eval(gen,%[]);

[A*X^2+B*X*Y+C*Y^2+D*X+E*Y+1 = 0, A+(3/2)*B+(9/4)*C+D+(3/2)*E+1 = 0, (1/4)*A-(1/2)*B+C-(1/2)*D+E+1 = 0, (16/9)*A-B+(9/16)*C+(4/3)*D-(3/4)*E+1 = 0, A+B+C-D-E+1 = 0]

[[A = -6*(1321*X*Y-3028*Y^2+471*X+873*Y+3051)/(18306*X^2+5401*X*Y-17332*Y^2-6054*X+12429*Y), B = (7926*X^2-2144*Y^2-3455*X+3836*Y-5401)/(18306*X^2+5401*X*Y-17332*Y^2-6054*X+12429*Y), C = -4*(4542*X^2-536*X*Y-2171*X+1844*Y-4333)/(18306*X^2+5401*X*Y-17332*Y^2-6054*X+12429*Y), D = (2826*X^2+3455*X*Y-8684*Y^2+3651*Y+6054)/(18306*X^2+5401*X*Y-17332*Y^2-6054*X+12429*Y), E = (5238*X^2-3836*X*Y+7376*Y^2-3651*X-12429)/(18306*X^2+5401*X*Y-17332*Y^2-6054*X+12429*Y)]]

-6*(1321*X*Y-3028*Y^2+471*X+873*Y+3051)*x^2/(18306*X^2+5401*X*Y-17332*Y^2-6054*X+12429*Y)+(7926*X^2-2144*Y^2-3455*X+3836*Y-5401)*x*y/(18306*X^2+5401*X*Y-17332*Y^2-6054*X+12429*Y)-4*(4542*X^2-536*X*Y-2171*X+1844*Y-4333)*y^2/(18306*X^2+5401*X*Y-17332*Y^2-6054*X+12429*Y)+(2826*X^2+3455*X*Y-8684*Y^2+3651*Y+6054)*x/(18306*X^2+5401*X*Y-17332*Y^2-6054*X+12429*Y)+(5238*X^2-3836*X*Y+7376*Y^2-3651*X-12429)*y/(18306*X^2+5401*X*Y-17332*Y^2-6054*X+12429*Y)+1 = 0

Check it out. Conic to pass through (-1/2,-3/4).

> 

con1:=eval(con,{X=-1/2,Y=-3/4});
c1:=implicitplot(con1,x=-2..2,y=-2..2,colour=black):display(c1,p1);

(66/109)*x^2+(200/327)*x*y-(512/327)*y^2-(5/109)*x+(76/109)*y+1 = 0

So we just need to systematically vary X,Y

> 

controlpts:=[seq([X,0.7*X],X=-2.1..2,0.5)];
c:=map(i->eval(con,{X=i[1],Y=i[2]}),controlpts): #list of conics

[[-2.1, -1.47], [-1.6, -1.12], [-1.1, -.77], [-.6, -.42], [-.1, -0.7e-1], [.4, .28], [.9, .63], [1.4, .98], [1.9, 1.33]]

> 

display(seq(implicitplot(c[i],x=-2..2,y=-2..2,colour=black),i=1..nops(c)),p1,axes=none);

> 

 

 

Download Pencil.mw

return as procedures...

dharr 8235
November 11 2020
1 0

I removed the with(orthopoly) since they aren't needed. I think you want to return the procedures themselves, so you need:

return psi,w; in each of functionA and functionB. Then it works.

Package.mw

packagecall.mw

 

pass Array...

dharr 8235
November 06 2020
0 0

As @acer pointed out in answer to your recent question here the way around this is to create the array outside the procedure and pass it as an argument to the procedure. I use compile quite a bit and usually end up passing many Arrays (or Vectors or Matrices); as @mmcdara says, there are lots of limitations.

convert pwlist...

dharr 8235
November 01 2020
1 0

convert(f,pwlist) puts the pieces in a list for easier manipulation.

show/hide...

dharr 8235
October 30 2020
1 8

Sounds like you are using worksheet mode rather than document mode. Try the view menu, choose show/hide contents and unclick input.

Complex Nonlinear Least Squares...

dharr 8235
October 08 2020
1 0

I wrote a Maple application for this, based on Maple's NLPSolve routine. It is available at the Maple Application Center

elementary...

dharr 8235
October 08 2020
1 2


 

> 

with(LinearAlgebra):with(ArrayTools):with(ListTools):
> 

E:=IdentityMatrix(5):
> 

Rowmult:=RowOperation(E,4,3);
Rowswap:=RowOperation(E,[1,3]);
Rowadd:=RowOperation(E,[2,3],5);

Rowmult := Matrix(5, 5, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1, (3, 4) = 0, (3, 5) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 3, (4, 5) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 1})

Rowswap := Matrix(5, 5, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 1, (1, 4) = 0, (1, 5) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (3, 1) = 1, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1, (4, 5) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 1})

Rowadd := Matrix(5, 5, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 5, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1, (3, 4) = 0, (3, 5) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1, (4, 5) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 1})

> 

Elementary:=proc(M::Matrix(square))
  local c,d,N;
  uses LinearAlgebra,ListTools,ArrayTools;
  N:=RowDimension(M);
  d:=Collect([seq(Diagonal(M))]);      #count things on diagonal
  if IsMatrixShape(M,diagonal)
     and nops(d)=2
     and member([1,N-1],d) then return true
  end if; #diagonal with N-1 1s and something else (row multiply)
  c:=Collect([seq(M)]);                #count zeroes and 1s and others
  if d=[[1,N]]
     and member([0,N^2-N-1],c) then return true
  end if; #1s on diagonal and one nonzero offdiagonal (row addition)   
  if nops(c)<>2 then return false end if; # row swap is a (0,1) matrix
  if not (c[1]=[0,N^2-N] and c[2]=[1,N]
       or c[2]=[0,N^2-N] and c[1]=[1,N])
     then return false
  end if; #not N ones and the rest zeroes
  if IsMatrixShape(M,symmetric)
     and Trace(M)=N-2
     and {seq(AddAlongDimension(M,1))}={1} #row sums 1
     then return true
   end if;                       # row swap matrix
  false  
end proc:
> 

Elementary(Rowmult);Elementary(Rowadd);Elementary(Rowswap);

true

true

true

> 

 

 

Download Elementary.mw

entries...

dharr 8235
October 07 2020
0 0

Elaborating on @nm's response, you could do something like this. (worksheet not displaying right now)

Download Numeric.mw

fraction-free...

dharr 8235
October 05 2020
1 3

There are fraction free options for GaussianElimination, LUDecopmposition and Determinant (maybe others), but not for Eigenvectors. You can normalize a vector of rationals to be fraction free with ff:=v->ilcm(seq(denom~(v)))*v;

uninteresting solution...

dharr 8235
September 22 2020
1 0

I corrected D^2 to D@@2, and knf = .8154646474. has two decimal points, so I deleted the second. Then it works, but your f and theta are just constant for the boundary conditions you have chosen.
 

> 

restart;
> 

with(plots):
> 

ODES := (diff(f(eta), `$`(eta, 4)))/((1-phi1)^2.5*(1-phi2)^2.5*((1-phi2)*(1-phi1+phi1*rhos1/rhosf)+phi2*rhos2/rhosf))+S*(f(eta)*(diff(f(eta), `$`(eta, 3)))-3*(diff(f(eta), `$`(eta, 2)))-eta*(diff(f(eta), `$`(eta, 3)))-(diff(f(eta), eta))*(diff(f(eta), `$`(eta, 2)))) = 0,

(khnf/kf+(4/3)*R)*(diff(theta(eta), `$`(eta, 2)))/((1-phi2)*(1-phi1+phi1*rhos1*cp1/(rhosf*cpf))+phi2*rhos2*cp2/(rhosf*cpf))+S*Pr*(f(eta)*(diff(theta(eta), eta))-eta*(diff(theta(eta), eta))-gamma*(eta^2*(diff(theta(eta), `$`(eta, 2)))-2*eta*f(eta)*(diff(theta(eta), `$`(eta, 2)))-eta*(diff(f(eta), eta))*(diff(theta(eta), eta))+f(eta)*(diff(f(eta), eta))*(diff(theta(eta), eta))+f(eta)^2*(diff(theta(eta), `$`(eta, 2))))) = 0;
 

(diff(diff(diff(diff(f(eta), eta), eta), eta), eta))/((1-phi1)^2.5*(1-phi2)^2.5*((1-phi2)*(1-phi1+phi1*rhos1/rhosf)+phi2*rhos2/rhosf))+S*(f(eta)*(diff(diff(diff(f(eta), eta), eta), eta))-3*(diff(diff(f(eta), eta), eta))-eta*(diff(diff(diff(f(eta), eta), eta), eta))-(diff(f(eta), eta))*(diff(diff(f(eta), eta), eta))) = 0, (khnf/kf+(4/3)*R)*(diff(diff(theta(eta), eta), eta))/((1-phi2)*(1-phi1+phi1*rhos1*cp1/(rhosf*cpf))+phi2*rhos2*cp2/(rhosf*cpf))+S*Pr*(f(eta)*(diff(theta(eta), eta))-eta*(diff(theta(eta), eta))-gamma*(eta^2*(diff(diff(theta(eta), eta), eta))-2*eta*f(eta)*(diff(diff(theta(eta), eta), eta))-eta*(diff(f(eta), eta))*(diff(theta(eta), eta))+f(eta)*(diff(f(eta), eta))*(diff(theta(eta), eta))+f(eta)^2*(diff(diff(theta(eta), eta), eta)))) = 0

> 

BCS:= f(0) = 0, ((D@@2)(f))(0) = 0, (D(theta))(0) = 0, f(1) = 0, (D(f))(1) = 0, theta(1) = 1;

 
 

f(0) = 0, ((D@@2)(f))(0) = 0, (D(theta))(0) = 0, f(1) = 0, (D(f))(1) = 0, theta(1) = 1

> 

params:={phi1 = .1, phi2 = .1, rhos1 = 2720, rhos2 = 2810, rhosf = 997.1, khnf = 1.083061737, kf = .613, cp1 = 893, cp2 = 960, cpf = 4179, Pr = 6.2, knf =.8154646474, S=0.5,R=0.5, gamma=0.5};

{Pr = 6.2, R = .5, S = .5, cp1 = 893, cp2 = 960, cpf = 4179, gamma = .5, kf = .613, khnf = 1.083061737, knf = .8154646474, phi1 = .1, phi2 = .1, rhos1 = 2720, rhos2 = 2810, rhosf = 997.1}

> 

sol:=dsolve(eval({ODES,BCS},params),numeric);

proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(8, {(1) = .0, (2) = .14285714285714277, (3) = .28571428571428553, (4) = .4285714285714284, (5) = .5714285714285713, (6) = .7142857142857142, (7) = .857142857142857, (8) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(8, 6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = 1.0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = 1.0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = 1.0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = 1.0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = 1.0, (5, 6) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = 1.0, (6, 6) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = 1.0, (7, 6) = .0, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0, (8, 5) = 1.0, (8, 6) = .0}, datatype = float[8], order = C_order); YP := Matrix(8, 6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = .0, (7, 6) = .0, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0, (8, 5) = .0, (8, 6) = .0}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(8, {(1) = .0, (2) = .14285714285714277, (3) = .28571428571428553, (4) = .4285714285714284, (5) = .5714285714285713, (6) = .7142857142857142, (7) = .857142857142857, (8) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(8, 6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = .0, (7, 6) = .0, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0, (8, 5) = .0, (8, 6) = .0}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[8] elif outpoint = "order" then return 2 elif outpoint = "error" then return HFloat(-0.0) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [6, 8, [f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta), diff(diff(diff(f(eta), eta), eta), eta), theta(eta), diff(theta(eta), eta)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[8] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[8] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(6, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(8, 6, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(6, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(8, 6, X, Y, outpoint, yout, L, V) end if; [eta = outpoint, seq('[f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta), diff(diff(diff(f(eta), eta), eta), eta), theta(eta), diff(theta(eta), eta)]'[i] = yout[i], i = 1 .. 6)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[8] elif outpoint = "order" then return 2 elif outpoint = "error" then return HFloat(-0.0) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [6, 8, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[8] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[8] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(8, 6, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(6, {(1) = 0., (2) = 0., (3) = 0., (4) = 0., (5) = 0., (6) = 0.}); `dsolve/numeric/hermite`(8, 6, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 6)] end proc, (2) = Array(0..0, {}), (3) = [eta, f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta), diff(diff(diff(f(eta), eta), eta), eta), theta(eta), diff(theta(eta), eta)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [eta = res[1], seq('[f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta), diff(diff(diff(f(eta), eta), eta), eta), theta(eta), diff(theta(eta), eta)]'[i] = res[i+1], i = 1 .. 6)] catch: error end try end proc

> 

odeplot(sol,[[eta,f(eta)],[eta,theta(eta)]],0..1,color=[red,blue]);

> 

 


 

Download dsolve.mw

first two terms...

dharr 8235
September 10 2020
0 0

series(HB1,r,3); gives the constant term and the term in r^2, but after that Maple seems to have a problem (same result with taylor)

adaptive=false...

dharr 8235
September 06 2020
0 1

The calculation is working, it is just taking forever. This can be seen by doing only 3 points, which produces a plot.

plot(evalf(eval(s__avg,[f=20e3,D_=0.2])),t=100e-6..200e-6,adaptive=false,sample=[100e-6,150e-6,200e-6],style=point,symbolsize=20);

Plotting with adaptive=false, numpoints=49 works, but you don't have good resolution on your steps. You can play around with the tradoff of accuracy and time.

You are numerically approximating an integral many times. It seems like normally you would use a sampled signal, and then the sliding average would reduce to a sum (add in Maple) and would be a much simpler calculation.

sliding.mw

display...

dharr 8235
September 02 2020
1 0

display combines plots. (For future reference it is helpful if you directly upload your worksheet using the big green up arrow.)
 

> 

restart;
> 

with(plots):
> 

alpha := 0.1;
beta := 0.1;
mu := 0.5;
u := 0.5;
v := 1;
gamma = 0.1;
 

.1

.1

.5

.5

1

gamma = .1

> 

Eq := -2*sigma*alpha*beta*mu*u - sigma*alpha*beta^2*u*v - 2.0*mu*alpha*beta^2*u*v - gamma^2*k^2*alpha*beta*u + beta*sigma^2*v + sigma^3 + 2*mu*sigma^2 - alpha*beta*sigma^2*u + sigma*beta^2*u*v + gamma^2*k^2*sigma + 2.0*mu*beta*sigma*v + 2.0*mu*beta^2*u*v;
p1:=implicitplot(Eq, k = 0 .. 10, sigma = -0.1 .. 0.1):

0.995e-1*sigma+0.4500e-2-0.5e-2*gamma^2*k^2+1.095*sigma^2+sigma^3+gamma^2*k^2*sigma

[k1,0], [k2,0], [k3,0], ...) , where k_{n}=2*n*Pi/L and the domain [0,L]

> 

L:=10.0;
p2:=pointplot([seq([2*n*Pi/L,0],n=0..15)],color=red,symbolsize=15):
display(p1,p2);

10.0

> 

 


 

Download plots.mw

x(t), y(t) in DEs...

dharr 8235
August 27 2020
1 8

In the DEs x and y should be x(t) and y(t) (I thought this was always true in Maple):

sys:=D(x)(t)=v1(x(t),y(t)),D(y)(t)=v2(x(t),y(t));

Note that the ditto used to be " but is now %, so that has to be changed also.

First 53 54 55 56 57 58 59 Last Page 55 of 81