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

identify...

dharr 8235
June 17 2019
1 0

identify(1.77245); gives sqrt(Pi)

asympt of RootOf...

dharr 8235
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 8235
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 8235
June 08 2019
2 6

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

combstruct...

dharr 8235
May 21 2019
1 1

 

> 

restart;
> 

with(combstruct):
> 

allstructs(Partition(10),size=2);

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

> 

 


 

Download partitions.mw

hard way...

dharr 8235
May 19 2019
0 0

[Edit: for S, the line

S:= (sqrt(2)/2^k)2*Matrix(M,M, (i,j)-> `if`(`and`(i::odd,j=1) ,-1/(2*i*(i-2)),0)):   

needs the  2 bolded above. Also, the loop you create S in can be discarded, you just need to set M and then the above line does it all]

For C, it's only partly a band matrix, and with many special cases in the first 2x2 block, I'd just make C this way

> 

restart;with(LinearAlgebra):
> 

M:=6;

6

Following needs to be multiplied by (1/2^k) to give C

> 

CC:=Matrix(M,M,(i,j)->if j=1 and i=1 then 1/2
                   elif i=1 and j=2 then 1/(2*sqrt(2))
                   elif j=1 and i=2 then -1/(8*sqrt(2))
                   elif j=1 and i>2 then -1/(2*sqrt(2)*i*(i-2))
                   elif j=i-1 then -1/(4*(i-2))
                   elif j=i+1 then -1/(4*i)
                   else 0
                   end if);

_rtable[18446744342674867014]

> 

 


 

Download MatrixC.mw

fracdiff...

dharr 8235
May 17 2019
0 3

For fractional differentiation in Maple use fracdiff. fracdiff(u(x,t),t,3/2) returns an integral. intsolve solves integral equations but complains this one is not linear, so I'm not sure where to go from here. Perhaps you know some methods for solving these eqns (laplace transforms perhaps?) that you could implement step by step.

math answer...

dharr 8235
May 14 2019
1 0

Premultiply by the 4x4 matrix with 1's on its reverse diagonal.
 

Download Matrix.mw

First 62 63 64 65 66 67 68 Last Page 64 of 81