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

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.

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

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.

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(1.77245); gives sqrt(Pi)

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

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

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


 

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

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

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

 

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 a symbolic solution, you can just use the rhs of the output of dsolve. And fprint should be fprintf.

wave.mw

[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

 

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