Advanced Search

dharr

Dr. David Harrington

838 Reputation

13 Badges

14 years, 323 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 professor of chemistry at the University of Victoria, BC, Canada, where my research areas are electrochemistry and surface science. I have been a user of Maple since about 1992.

MaplePrimes Activity


These are answers submitted by dharr

DFT...

dharr 838
July 11 2019
1 2

Fourier series usually means the infinite series, yours looks more like a discrete fourier transform. Search DFT - there is a version in the SignalProcessing package called DFT (or FFT for efficiency) and another in the DiscreteTransforms package called FourierTransform. The formulations are not exactly as you have, but run from 0 to N-1, and have the exponential negative for the forward transform and positive for the inverse transform.

zero...

dharr 838
July 10 2019
3 2

If I leave tau undefined and change "evalf" to "simplify" so there are no floats, the integral evaluates to zero. Is there a reason to thnk this value is wrong?

070919_ThetaIntegral.mw

f1,1 and f11...

dharr 838
July 10 2019
1 0

The f_1,1 in the first line is diiferent from the f_11 later on. If you make them all f_11 it works.

2-D notation problem?...

dharr 838
July 09 2019
1 0

Since you didn't upload your worksheet, it is hard to diagnose, but it works here in 1-D:
 

> 

restart;
> 

with(PDEtools):
> 

declare(u(x,t));
U:=diff_table(u(x,t));

u(x, t)*`will now be displayed as`*u

table( [( ) = u(x, t) ] )

> 

PDE := U[x,x]-U[t];

diff(diff(u(x, t), x), x)-(diff(u(x, t), t))

> 

Infinitesimals(PDE);

[_xi[x](x, t, u) = 0, _xi[t](x, t, u) = 1, _eta[u](x, t, u) = 0], [_xi[x](x, t, u) = 1, _xi[t](x, t, u) = 0, _eta[u](x, t, u) = 0], [_xi[x](x, t, u) = 0, _xi[t](x, t, u) = 0, _eta[u](x, t, u) = u], [_xi[x](x, t, u) = (1/2)*x, _xi[t](x, t, u) = t, _eta[u](x, t, u) = 0], [_xi[x](x, t, u) = -2*t, _xi[t](x, t, u) = 0, _eta[u](x, t, u) = x*u], [_xi[x](x, t, u) = 0, _xi[t](x, t, u) = 0, _eta[u](x, t, u) = exp(_c[1]*t)*exp(_c[1]^(1/2)*x)], [_xi[x](x, t, u) = 0, _xi[t](x, t, u) = 0, _eta[u](x, t, u) = exp(_c[1]*t)/exp(_c[1]^(1/2)*x)], [_xi[x](x, t, u) = (1/2)*x*t, _xi[t](x, t, u) = (1/2)*t^2, _eta[u](x, t, u) = -(1/8)*x^2*u-(1/4)*u*t]

> 

 


 

Download Infinitesimals.mw

identify...

dharr 838
June 17 2019
1 0

identify(1.77245); gives sqrt(Pi)

asympt of RootOf...

dharr 838
June 16 2019
2 5
> 

restart;
> 

eq:=((2*n-1)*x + 2*n + 1)*x^n = 1 - x;

((2*n-1)*x+2*n+1)*x^n = 1-x

> 

solve(eq,x);
asympt(%,n);
ass:=convert(%,polynom);

RootOf(2*_Z^n*n*_Z+2*_Z^n*n-_Z^n*_Z+_Z^n+_Z-1)

1-LambertW(4*n^2)/n+(1/2)*LambertW(4*n^2)^2/n^2-(1/12)*LambertW(4*n^2)^2*(2*LambertW(4*n^2)^2+3*LambertW(4*n^2)+3)/((LambertW(4*n^2)+1)*n^3)+(1/24)*LambertW(4*n^2)^3*(LambertW(4*n^2)^2+3*LambertW(4*n^2)+6)/((LambertW(4*n^2)+1)*n^4)-(1/1440)*LambertW(4*n^2)^3*(12*LambertW(4*n^2)^5+89*LambertW(4*n^2)^4+312*LambertW(4*n^2)^3+435*LambertW(4*n^2)^2+270*LambertW(4*n^2)+90)/((LambertW(4*n^2)+1)^3*n^5)+O(1/n^6)

1-LambertW(4*n^2)/n+(1/2)*LambertW(4*n^2)^2/n^2-(1/12)*LambertW(4*n^2)^2*(2*LambertW(4*n^2)^2+3*LambertW(4*n^2)+3)/((LambertW(4*n^2)+1)*n^3)+(1/24)*LambertW(4*n^2)^3*(LambertW(4*n^2)^2+3*LambertW(4*n^2)+6)/((LambertW(4*n^2)+1)*n^4)-(1/1440)*LambertW(4*n^2)^3*(12*LambertW(4*n^2)^5+89*LambertW(4*n^2)^4+312*LambertW(4*n^2)^3+435*LambertW(4*n^2)^2+270*LambertW(4*n^2)+90)/((LambertW(4*n^2)+1)^3*n^5)

> 

n:=1000;
fsolve(eq,x=0..1);
evalf(ass);

1000

.9874167130

.9874167131

> 

 


 

Download limit.mw

ad hoc solution...

dharr 838
June 15 2019
1 1

If I understand correctly,you know B,and you only know the eigenvalues relative to the most negative one. If so, then you can work out V3 and V6 from the symbolic solutions for the eigenvalues - I'm also assuming you know which are degenerate. Your matrix has a lot of structure, so it might be possible to generalize this.

Download HinderedRotor2.mw

missing multiplication and missing term...

dharr 838
June 08 2019
1 5

You do not have a space between n and (n+1). Insert a space or multiplication sign between them and you will get a different solution, presumable the one you want. However, (see ?LegendreP), your differential equation needs another term to qualify - adding that leads to a solution in terms of Legendre functions.


 

NULL

NULL

> 

restart
> 

ode := (-x^2+1)*(diff(y(x), x, x))+n*(n+1)*y(x) = 0

(-x^2+1)*(diff(diff(y(x), x), x))+n*(n+1)*y(x) = 0

> 

sol[1] := rhs(dsolve(ode))

_C1*(-x^2+1)*hypergeom([1+(1/2)*n, 1/2-(1/2)*n], [1/2], x^2)+_C2*(-x^3+x)*hypergeom([1-(1/2)*n, 3/2+(1/2)*n], [3/2], x^2)

> 

`assuming`([convert(sol[1], StandardFunctions)], [n::posint])

_C1*(-x^2+1)*hypergeom([1+(1/2)*n, 1/2-(1/2)*n], [1/2], x^2)+_C2*(-x^3+x)*hypergeom([1-(1/2)*n, 3/2+(1/2)*n], [3/2], x^2)

> 

ode2 := (-x^2+1)*(diff(y(x), x, x))-2*x*(diff(y(x), x))+n*(n+1)*y(x) = 0

(-x^2+1)*(diff(diff(y(x), x), x))-2*x*(diff(y(x), x))+n*(n+1)*y(x) = 0

> 

sol[2] := rhs(dsolve(ode2))

_C1*LegendreP(n, x)+_C2*LegendreQ(n, x)

NULL


 

Download convert-Legendre.mw

multiple roots...

dharr 838
May 27 2019
0 0


 

> 

restart;
> 

eqs := [a+b+c=10, a/b=b/c, a*b=6];

[a+b+c = 10, a/b = b/c, a*b = 6]

> 

fsolve(eqs); ## finds one real root

{a = 1.935562277, b = 3.099874424, c = 4.964563299}

> 

solve(eqs);  #has quartic, so potentially 4 roots

{a = -(1/6)*RootOf(_Z^4+6*_Z^2-60*_Z+36)^3-RootOf(_Z^4+6*_Z^2-60*_Z+36)+10, b = RootOf(_Z^4+6*_Z^2-60*_Z+36), c = (1/6)*RootOf(_Z^4+6*_Z^2-60*_Z+36)^3}

> 

evalf(allvalues(solve(eqs)));

{a = 9.311003375, b = .644398865, c = 0.4459776043e-1}, {a = 1.935562295, b = 3.099874421, c = 4.964563284}, {a = -.62328281-1.268492836*I, b = -1.872136643+3.810135335*I, c = 12.49541945-2.541642499*I}, {a = -.62328281+1.268492836*I, b = -1.872136643-3.810135335*I, c = 12.49541945+2.541642499*I}

> 

 


 

Download solve.mw

dsolve with parameters...

dharr 838
May 27 2019
1 1

You need to tell dsolve about the parameters, and then declare their values before odeplot. This can get you separate odeplots (eg., vs time) for different parameters, see attached. You could then combine them with display to get a series of plots at different parameter values, as below.

[Had problem with upload - removed output]

Briefly: don't have a value for E at the time you use dsolve, but add parameters=[E] to dsolve.

Then set a value before you odeplot by

Q(parameters=[2.5])

NULL

> 

restart; with(linalg); with(VectorCalculus); with(DEtools); with(plots); r := .9; K := 10; beta := .5; g := 0.3e-1; alpha[1] := .4; alpha[2] := 2.2; alpha[3] := 4; d[1] := .1; d[2] := .1; q := 1; s := .46; delta := 1.5; b := 1.2; p := 1; c := 0.1e-1; mu := 0.25e-1; P0 := 4; M0 := 5; N0 := 3; L0 := 2; T := 20

Model 1;

 
> 

dP := r*P(t)*(1-P(t)/K)+g*P(t)-beta*P(t); dM := beta*P(t)-q*E*M(t)-d[1]*M(t); dN := -s*N(t)+delta*N(t)*M(t)/(M(t)+b*N(t)); dL := (alpha[1]-alpha[2]*L(t)/(alpha[3]+M(t)))*L(t)-d[2]*L(t)
> 

satu := diff(P(t), t) = dP; dua := diff(M(t), t) = dM; tiga := diff(N(t), t) = dN; empat := diff(L(t), t) = dL
> 

pdb := satu, dua, tiga, empat; fcns := {L(t), M(t), N(t), P(t)}

> 

Q := dsolve({pdb, L(0) = L0, M(0) = M0, N(0) = N0, P(0) = P0}, fcns, type = numeric, method = rkf45, maxfun = 500000, parameters = [E])
> 

Q(parameters = [2.5])

> 

p25 := odeplot(Q, [[t, P(t), color = blue], [t, M(t), color = green], [t, N(t), color = red], [t, L(t), color = gold]], t = 0 .. T, numpoints = 100000, thickness = 2); display(p25)
> 

Q(parameters = [4.0]); p40 := odeplot(Q, [[t, P(t), color = blue], [t, M(t), color = green], [t, N(t), color = red], [t, L(t), color = gold]], t = 0 .. T, numpoints = 100000, thickness = 2); display(p40)

Combine plots

> 

display(p25, p40)

NULL

NULL


 

Download captive_breeding.mw

solve tomleslie's recurrence...

dharr 838
May 26 2019
1 0

Or, following on from @tomleslie's recurrence:
 

> 

  restart;
  with(genfunc):
  p:=[2,3,4,6,8,9,10,12,14,15,16,18];
#
# Produce recurrence relation
# and solve it
#
  rgf_findrecur(6, p, f, j);
  rsolve({%,f(1)=2,f(2)=3,f(3)=4,f(4)=6},f(n));

[2, 3, 4, 6, 8, 9, 10, 12, 14, 15, 16, 18]

f(j) = 2*f(j-1)-2*f(j-2)+2*f(j-3)-f(j-4)

-((1/4)*I)*I^n+((1/4)*I)*(-I)^n+(3/2)*n

> 

 


 

Download sequence2.mw

 

sequence...

dharr 838
May 26 2019
1 0
> 

restart;
> 

s:=[2,3,4,6,8,9,10,12,14,15,16,18];
i:=[$1..nops(s)];

[2, 3, 4, 6, 8, 9, 10, 12, 14, 15, 16, 18]

[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]

Note

> 

seq(s[n],n=2..12,2); #multiples of 3
seq(s[n],n=1..12,4); #increments of 6
seq(s[n],n=3..12,4); #increments of 6

3, 6, 9, 12, 15, 18

2, 8, 14

4, 10, 16

> 

f:=n->if n::even then 3*n/2
      elif n mod 4 = 1 then 3*(n-1)/2+2
      else 3*(n-1)/2+1 end if;
map(f,i);

proc (n) options operator, arrow; if n::even then (3/2)*n elif `mod`(n, 4) = 1 then (3/2)*n+1/2 else (3/2)*n-1/2 end if end proc

[2, 3, 4, 6, 8, 9, 10, 12, 14, 15, 16, 18]

Or, at least for the given numbers

> 

g:=CurveFitting[PolynomialInterpolation](i,s,n);

(1/2494800)*n^11-(1/37800)*n^10+(17/22680)*n^9-(1/84)*n^8+(2927/25200)*n^7-(431/600)*n^6+(31901/11340)*n^5-(2579/378)*n^4+(140776/14175)*n^3-(13304/1575)*n^2+(35611/6930)*n

> 

eval(g,n=10);

15

> 

 


 

Download sequence.mw

from symbolic solution...

dharr 838
May 25 2019
0 1

From a symbolic solution, you can just use the rhs of the output of dsolve. And fprint should be fprintf.

wave.mw

output Array...

dharr 838
May 24 2019
1 2

[Edit: this was posted before OP uploaded worksheet]

I'm not sure if the output=Array option was introduced before Maple 13, but if it was, this should work.
 

> 

restart;
> 

edo:=diff(u(t),t,t)=-3*u(t);
ics:=D(u)(0)=2,u(0)=0;

diff(diff(u(t), t), t) = -3*u(t)

(D(u))(0) = 2, u(0) = 0

Times you want the output at

> 

times:=Array([seq(0.1*i,i=0..20)]);

times := Vector(4, {(1) = ` 1 .. 21 `*Array, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

> 

ans:=dsolve({edo,ics},u(t),numeric,output=times);

_rtable[18446744576662716894]

Write t,u(t) pairs to a file

> 

writedata(cat(currentdir(),"/testfile.txt"),convert(ans[2,1][..,1..2],listlist));
> 

 


 

Download edo.mw

 

function of x and y...

dharr 838
May 22 2019
0 0

plot3d of Re(Q) can plot Re(Q) as a function of two variables, say x and y, so the t range should be omitted and t and a given values (note gamma has a defined value of 0.577215 see ?gamma - if you didn't intend that then choose a different name). Then you can make a plot as below

Download Q2.mw

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