Advanced Search

dharr

Dr. David Harrington

8205 Reputation

22 Badges

20 years, 336 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

multiple roots...

dharr 8205
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 8205
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 8205
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 8205
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 8205
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 8205
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 8205
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 8205
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 8205
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 8205
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 8205
May 14 2019
1 0

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

Download Matrix.mw

diff...

dharr 8205
May 08 2019
2 1

The 1 or 2 means differential with respect to the first or second variable, so [2,2] means with respect to the (unnamed) second variable and then again with respect to that variable. The (x,0) means that the second variable is then set to zero. 


 

> 

a:=(D[2, 2](T))(x, 0) = 0; b:= (D[1, 1](T))(x, 0) = 0;

(D[2, 2](T))(x, 0) = 0

(D[1, 1](T))(x, 0) = 0

> 

convert(a,Diff);convert(b,Diff);

eval(Diff(T(x, t1), t1, t1), {t1 = 0}) = 0

Diff(T(x, 0), x, x) = 0

> 

 


 

Download diff.mw

deriving the recurrence...

dharr 8205
May 07 2019
2 3

Let e(n) be sequences with n numbers ending in an even number, o(n) be those with n numbers ending in an odd number.

We start with e(1)=2 (0 or 2) and o(n)=1 (1). We can add a 1 to the end of an even or and odd sequence, o(n+1)=e(n)+o(n); we can only add 0 or 2 to the end of an odd sequence, so e(n+1)=2*o(n)

> 

ans:=rsolve({o(n+1)=e(n)+o(n),e(n+1)=2*o(n),e(1)=2,o(1)=1},{e(n),o(n)});

{e(n) = -(2/3)*(-1)^n+(2/3)*2^n, o(n) = (1/3)*(-1)^n+(2/3)*2^n}

> 

total:=unapply(eval(e(n)+o(n),ans),n);

proc (n) options operator, arrow; -(1/3)*(-1)^n+(4/3)*2^n end proc

> 

map(total,[$(1..6)]);

[3, 5, 11, 21, 43, 85]

> 

 


 

Download ternary.mw

intersection is a line...

dharr 8205
May 06 2019
0 1

The variable in LinearSolve means there is a family of solutions, a line, which can be plotted with spacecurve. The planes can be plotted in various ways, e.g. with implicitplot3d

problem.mw

table width...

dharr 8205
May 06 2019
1 2

Not sure if I understand what you want. If you right-click and find table properties, you can change the width of the table to be say 50% of the page width. This seems to be respected when you export to .pdf

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