Advanced Search

dharr

Dr. David Harrington

8200 Reputation

22 Badges

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

not autonomous...

dharr 8200
August 08 2024
3 1

"A system of two first order differential equations produces a direction field plot, provided the system is determined to be autonomous. In addition, a single first order differential equation produces a direction field (as it can always be mapped to a system of two first order autonomous differential equations). A system is determined to be autonomous when all terms and factors, other than the differential, are free of the independent variable. "

no solutions in range...

dharr 8200
August 08 2024
1 0

The syntax is correct, just there are no solutions in this range. Works for -2..2

> 

restart
> 

f := sin(x[1]+x[2])-exp(x[1])*x[2]; g := x[1]^2-x[2]-2; cond := {seq(x[i] = -1 .. 1, i = 1 .. 2)}; fsolve({f, g}, cond)

sin(x[1]+x[2])-exp(x[1])*x[2]

x[1]^2-x[2]-2

{x[1] = -1 .. 1, x[2] = -1 .. 1}

fsolve({sin(x[1]+x[2])-exp(x[1])*x[2], x[1]^2-x[2]-2}, {x[1], x[2]}, {x[1] = -1 .. 1, x[2] = -1 .. 1})

> 

plot3d([f, g], x[1] = -2 .. 2, x[2] = -2 .. 2)

> 

cond2 := {seq(x[i] = -2 .. 2, i = 1 .. 2)}; fsolve({f, g}, cond2)

{x[1] = -2 .. 2, x[2] = -2 .. 2}

{x[1] = -.6687012050, x[2] = -1.552838698}

> 

NULL

Download fsolve.mw

One way...

dharr 8200
August 08 2024
1 0

Seems like overkill, but here is one way. As for why, ...
(Note it is complexfreqvar that was s, not statevariable.)

> 

restart;
> 

foo:=proc()
 local stemp;
 if assigned(':-s') then
   stemp := :-s;
   :-s := ':-s';
   DynamicSystems:-SystemOptions('complexfreqvar'=':-sv');
   :-s := stemp;
 else
    DynamicSystems:-SystemOptions('complexfreqvar'=':-sv');
 end if;
 :-sv
end proc:
> 

s:="test";
foo();
s;

"test"

sv

"test"

> 

 

Download dynamic.mw

solve for z...

dharr 8200
August 08 2024
1 1 
sol := solve({eq1, eq2, eq3, x + y + z = 1}, {a, b, c, z}, explicit);

faster...

dharr 8200
August 08 2024
1 2

I have f(0, 12) taking 44 s, so you might get a result for f(0,30) in managable time.

For future reference please upload your worksheet using the green up-arrow in the Mapleprimes editor.

> 

restart

Use rationals

> 

n1 := 423; x := 16; n2 := 81; y := 35; s1 := 1/10; b1 := 7; beta1 := 21/10; s2 := 1/10; b2 := 7; beta2 := 21/10;

423

16

81

35

1/10

7

21/10

1/10

7

21/10

> 

A := (h, v, r) -> add(add(binomial(n1 - x, j)*binomial(n2 - y, l)*(-1)^(j + l)*GAMMA(s1 + h + 1)*GAMMA(s2 + v + 1)*(-1)^(h + v)*(beta1 - 1)^h*(beta2 - 1)^v/(h!*GAMMA(s1 + 1)*v!*GAMMA(s2 + 1)*(x + b1*s1 + j + b1*h + y + b2*s2 + l + b2*v)*(r + x + b1*s1 + j + b1*v)), j = 0 .. n1 - x), l = 0 .. n2 - y);

Warning, (in A) `l` is implicitly declared local

Warning, (in A) `j` is implicitly declared local

proc (h, v, r) local l, j; options operator, arrow; add(add(binomial(n1-x, j)*binomial(n2-y, l)*(-1)^(l+j)*GAMMA(s1+h+1)*GAMMA(s2+v+1)*(-1)^(h+v)*(beta1-1)^h*(beta2-1)^v/(factorial(h)*GAMMA(s1+1)*factorial(v)*GAMMA(s2+1)*(x+b1*s1+j+b1*h+y+b2*s2+l+b2*v)*(r+x+b1*s1+j+b1*v)), j = 0 .. n1-x), l = 0 .. n2-y) end proc

Result is a rational, but very small value, so will have numerical issues if you use floats

> 

A(2,2,0);evalf(%);

11727199469241514779864393729287999595266717038898384405454822453889989726964576758757729551099323484694782170822448350244634045448648056150840927233527350118380124119264858615283201523602923810903830006171445640638533014964609388484372553536841310967853824558690932280088151544584971199415788181111551306532754880119508566210291084045598562471540440317095151211911113829813243986384020095260781211101363284506805774438303489939264952762167438202379620004245323297602579179457497328410039101743209644168235298998961095392278935981336337925654946855818892100272459692432824817768478928223155226897063053135874218725475831316965225841532156062303907670391646481346412180543470550738955412983380825233032396184376956505282354537912590976013025195012646732190615857238658037431378033943474292755126953125/6382265268317256219633710058998499818340873590563036702803822950484199699636302738462052037083809489945501654638778705432982835407262106422716055822223400225220629848888603203921980303115472012139298206087469920885254016501522878034151924440738350864243996392880003254552368644323922787997879524243559004914030851003740366052093558872264913393738584398643237741388376227019971696955206515106578177861348809414199987981935214749218771677333981945489142846776730080142869549055976924435701756372241817375336920507491648219754129346141714754707551648893025832772408084923821700316568058764587634659823042662378587917652228685077432247136823506786967161644710468316981814899035776299987989971099597861998466613316946045883428522261542320902882571202585668513207901611181230727412654811333329187224660577277205763279177427978308304241574151519438327049907892296647014672715620294656

0.1837466632e-77

Compare with a floating point calculation!

> 

A(2, 2, 0.)

-0.3955624180e121

> 

B := (h, v, r) -> add(add(binomial(n1 - x, j)*binomial(n2 - y, l)*(-1)^(j + l)*GAMMA(s1 + v + 1)*GAMMA(s2 + h + 1)*(-1)^(h + v)*(beta1 - 1)^v*(beta2 - 1)^h/(v!*GAMMA(s1 + 1)*h!*GAMMA(s2 + 1)*(y + b2*s2 + l + b2*h + x + b1*s1 + j + b1*v)*(y + b2*s2 + l + b2*h - r)), j = 0 .. n1 - x), l = 0 .. n2 - y);

Warning, (in B) `l` is implicitly declared local

Warning, (in B) `j` is implicitly declared local

proc (h, v, r) local l, j; options operator, arrow; add(add(binomial(n1-x, j)*binomial(n2-y, l)*(-1)^(l+j)*GAMMA(s1+v+1)*GAMMA(s2+h+1)*(-1)^(h+v)*(beta1-1)^v*(beta2-1)^h/(factorial(v)*GAMMA(s1+1)*factorial(h)*GAMMA(s2+1)*(y+b2*s2+l+b2*h+x+b1*s1+j+b1*v)*(y+b2*s2+l+b2*h-r)), j = 0 .. n1-x), l = 0 .. n2-y) end proc

> 

f := (r, upto) -> sum(sum(A(h, v, r) + B(h, v, r), v = 0 .. upto), h = 0 .. upto);

proc (r, upto) options operator, arrow; sum(sum(A(h, v, r)+B(h, v, r), v = 0 .. upto), h = 0 .. upto) end proc

> 

s := CodeTools:-Usage(f(0,12)):

memory used=14.57GiB, alloc change=214.03MiB, cpu time=32.14s, real time=44.15s, gc time=24.25s

> 

s;

114249022902597666588491400529505235216742049616231180869631619194167698464302221521109360278934161699249264126826236228253701293270492644550135177145016349814982113455183820820735876815796463679992519208283503642450130209375877055834771405738601458875577451563602032933837943025767731232207854934759236939341322568521090401156938484109405395425106996068885180201446511126644125409694777192264543924174126949778005106553840262917702653247531868774015140088404720078587672388076732410686714770939922375245824000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000/4478112673625161634402912271274057087521425782459964471363322850631457597742398047074088546911131451856756368982860023346654177624065136219112870362974132591633772922976976016146648645581986021657025573763585023497503972196876316120665545360282282680322522472194461718178716824852465782992306776138047404097571796261065850414338472933710797700189622470162561730240651195333635833352277583626670180879730361251512986788399320104862652482661166302364468264845833967133426250993316566352774199352462247496894848151055424852777184106477710417868463737767180212637570654326336497085984319988974445914462490654328034321182352155451367268612522705340661230039301539950445376032066380305314154637436491880112866704742279027028488977660000585169277423676424285496535951865146024850581707926633898571343437487693846692653268159361127272552016275556640003360341921899571661034103268149587939236207824679329125781964757872174171342424591595209864084608515982176191174973912465653082613424079377178453298484585210385316658608895204147668705453968995722653400843313999068483592478297363991

> 

evalf(s);

0.2551276201e-55

> 

NULL

Download sum.mw

agreed...

dharr 8200
August 07 2024
2 1

Maple likes the expanded form and one has to "fake it" in order to see exactly what you want. These sorts of manipulations can usually be done, but are not easy and IMO not worth the effort. Maple likes two fewer parentheses - isn't that good?

> 

restart
> 

p := 9*csc(theta)^2+9

9*csc(theta)^2+9

Use of %* instead of * stops the automatic multiplication

> 

q := content(p)*%*primpart(p); value(q)

`%*`(9, csc(theta)^2+1)

9*csc(theta)^2+9

Less ugly but more complicated

> 

InertForm:-Display(content(p)*%*primpart(p), inert = false)

0, "%1 is not a command in the %2 package", _Hold, Typesetting

> 

NULL

Download content.mw

indexing error...

dharr 8200
August 07 2024
0 0

Your code uses u in three different ways: as a procedure, as indexed variables, and as a 2-dimensional Array(1..40, 0..4). These should be using different symbols. By the time you enter eval_derivatives, u is the Array, and the other uses are gone. Then you use u[0], which means the zeroth row of the Array. But the rows start at 1, so you get the error message.

Also, I see later in eval_derivatives, you use T1i[i, k - j], when T1i is a 1D Array.

generic...

dharr 8200
August 06 2024
2 0

By default, Maple assumes generic values of the coefficients. For special values, one needs to work harder. Here is one way to do it.

Edit: solve([eq1, eq2], [A__1, A__2, w], explicit) works also.

SolveTools:-PolynomialSystem([eq1, eq2], [A__1, A__2, w], explicit)
> 

restart
> 

eq1 := A__1*(w^2-3*w__0^2)+A__2*w__0^2;
eq2 := A__1*w__0^2+A__2*(w^2-w__0^2)NULL

A__1*(w^2-3*w__0^2)+A__2*w__0^2

A__1*w__0^2+A__2*(w^2-w__0^2)

If we treat it as a polynomial system

> 

SolveTools:-PolynomialSystem([eq1, eq2], [A__1, A__2, w], explicit)

{A__1 = 0, A__2 = 0, w = w}, {A__1 = -A__2*(1+2^(1/2)), A__2 = A__2, w = (2+2^(1/2))^(1/2)*w__0}, {A__1 = -A__2*(1+2^(1/2)), A__2 = A__2, w = -(2+2^(1/2))^(1/2)*w__0}, {A__1 = -A__2*(1-2^(1/2)), A__2 = A__2, w = (2-2^(1/2))^(1/2)*w__0}, {A__1 = -A__2*(1-2^(1/2)), A__2 = A__2, w = -(2-2^(1/2))^(1/2)*w__0}

> 

NULL

Download solve.mw

table...

dharr 8200
August 06 2024
2 1

I rewrote this as a table. You don't say if you want all possible p and q sets for each expression, or just one. The following shows both possibilities. As I told you before, you were not properly detecting duplicates because the same expression was appearing in slightly different forms. evalc is enough to make a canonical form here; simplify is more obvious but too expensive.

There are various changes that could be make depending on how you want the output.

> 

restart;
> 

result:=table(sparse=NULL): #sparse=NULL only required if accumulating
phi := (p__1, p__2, p__3, p__4, q__1, q__2, q__3, q__4, xi) -> evalc((p__1*cosh(q__1*xi) - p__2*sinh(q__2*xi))/(p__3*cosh(q__3*xi) + p__4*sinh(q__4*xi)));
for p__1 in [1, -1, I, -I] do
    for p__2 in [1, -1, I, -I] do
        for p__3 in [1, -1, I, -I] do
            for p__4 in [ 0] do
                for q__1 in [ I, -I] do
                    for q__2 in [ I, -I] do
                        for q__3 in [1] do
                            for q__4 in [ I, -I] do
                                # next line overwrites if we get the same result
                                result[phi(p__1, p__2, p__3, p__4, q__1, q__2, q__3, q__4, xi)] := [p__1, p__2, p__3, p__4, q__1, q__2, q__3, q__4];
                                # OR, next line accumulates all possibilities
                                #result[phi(p__1, p__2, p__3, p__4, q__1, q__2, q__3, q__4, xi)] ,= [p__1, p__2, p__3, p__4, q__1, q__2, q__3, q__4];
                            od:
                       od:
                    od:
                od:
            od:
         od:
     od:
od:

proc (p__1, p__2, p__3, p__4, q__1, q__2, q__3, q__4, xi) options operator, arrow; evalc((p__1*cosh(q__1*xi)-p__2*sinh(q__2*xi))/(p__3*cosh(q__3*xi)+p__4*sinh(q__4*xi))) end proc

All the expressions

> 

exprns:=[indices(result,'nolist')];
print("Number of possibilities is ",numelems(exprns));

[(-cos(xi)+sin(xi))/cosh(xi), -sin(xi)/cosh(xi)-I*cos(xi)/cosh(xi), -cos(xi)/cosh(xi)+I*sin(xi)/cosh(xi), I*(-cos(xi)+sin(xi))/cosh(xi), cos(xi)/cosh(xi)-I*sin(xi)/cosh(xi), I*(cos(xi)-sin(xi))/cosh(xi), cos(xi)/cosh(xi)+I*sin(xi)/cosh(xi), (cos(xi)+sin(xi))/cosh(xi), sin(xi)/cosh(xi)+I*cos(xi)/cosh(xi), -cos(xi)/cosh(xi)-I*sin(xi)/cosh(xi), sin(xi)/cosh(xi)-I*cos(xi)/cosh(xi), (cos(xi)-sin(xi))/cosh(xi), -sin(xi)/cosh(xi)+I*cos(xi)/cosh(xi), (-cos(xi)-sin(xi))/cosh(xi), I*(-cos(xi)-sin(xi))/cosh(xi), I*(cos(xi)+sin(xi))/cosh(xi)]

"Number of possibilities is ", 16

Find the corresponding p and q values. For just one

> 

result[(-cos(xi) + sin(xi))*I/cosh(xi)];
result[exprns[3]];

[-I, -1, 1, 0, -I, I, 1, -I]

[-I, -I, I, 0, -I, I, 1, -I]

For all of them

> 

[entries(result,pairs)];

[(-cos(xi)+sin(xi))/cosh(xi) = [-I, -1, I, 0, -I, I, 1, -I], -sin(xi)/cosh(xi)-I*cos(xi)/cosh(xi) = [-I, -I, 1, 0, -I, I, 1, -I], -cos(xi)/cosh(xi)+I*sin(xi)/cosh(xi) = [-I, -I, I, 0, -I, I, 1, -I], I*(-cos(xi)+sin(xi))/cosh(xi) = [-I, -1, 1, 0, -I, I, 1, -I], cos(xi)/cosh(xi)-I*sin(xi)/cosh(xi) = [-I, -I, -I, 0, -I, I, 1, -I], I*(cos(xi)-sin(xi))/cosh(xi) = [-I, -1, -1, 0, -I, I, 1, -I], cos(xi)/cosh(xi)+I*sin(xi)/cosh(xi) = [-I, -I, -I, 0, -I, -I, 1, -I], (cos(xi)+sin(xi))/cosh(xi) = [-I, -1, -I, 0, -I, -I, 1, -I], sin(xi)/cosh(xi)+I*cos(xi)/cosh(xi) = [-I, -I, -1, 0, -I, I, 1, -I], -cos(xi)/cosh(xi)-I*sin(xi)/cosh(xi) = [-I, -I, I, 0, -I, -I, 1, -I], sin(xi)/cosh(xi)-I*cos(xi)/cosh(xi) = [-I, -I, 1, 0, -I, -I, 1, -I], (cos(xi)-sin(xi))/cosh(xi) = [-I, -1, -I, 0, -I, I, 1, -I], -sin(xi)/cosh(xi)+I*cos(xi)/cosh(xi) = [-I, -I, -1, 0, -I, -I, 1, -I], (-cos(xi)-sin(xi))/cosh(xi) = [-I, -1, I, 0, -I, -I, 1, -I], I*(-cos(xi)-sin(xi))/cosh(xi) = [-I, -1, 1, 0, -I, -I, 1, -I], I*(cos(xi)+sin(xi))/cosh(xi) = [-I, -1, -1, 0, -I, -I, 1, -I]]

> 

NULL

Download generate_solution3.mw

Mapleprimes...

dharr 8200
August 05 2024
2 1

Sometimes the Mapleprimes worksheet display is broken for a time, and the files display incorrectly as you saw. But sometimes there is just something slightly different about a file that leads to problems, as is the case here. I opened your file and had the same poblem uploading, then deleted the last two plots and then redid them. Now it seems identical but can upload. 

This project discusses predator-prey system, particularly the Lotka-Volterra equations,which model the interaction between two sprecies: prey and predators. Let's solve the Lotka-Volterra equations numerically and visualize the results.

NULL

NULL

alpha := 1.0; beta := .1; g := 1.5; delta := 0.75e-1; ode1 := diff(x(t), t) = alpha*x(t)-beta*x(t)*y(t); ode2 := diff(y(t), t) = delta*x(t)*y(t)-g*y(t); eq1 := -beta*x*y+alpha*x = 0; eq2 := delta*x*y-g*y = 0; equilibria := solve({eq1, eq2}, {x, y}); print("Equilibrium Points: ", equilibria); initial_conditions := x(0) = 40, y(0) = 9; sol := dsolve({ode1, ode2, initial_conditions}, {x(t), y(t)}, numeric); eq_points := [seq([rhs(eq[1]), rhs(eq[2])], `in`(eq, equilibria))]

[[0., 0.], [20., 10.]]

plots[odeplot](sol, [[t, x(t)], [t, y(t)]], t = 0 .. 100, legend = ["Rabbits", "Wolves"], title = "Prey-Predator Dynamics", labels = ["Time", "Population"])

NULL

NULL

NULL

sol_plot := plots:-odeplot(sol, [[x(t), y(t)]], 0 .. 100, color = "blue"); equilibrium_plot := plots:-pointplot(eq_points, color = "red", symbol = solidcircle, symbolsize = 15); plots:-display([sol_plot, equilibrium_plot], title = "Phase Portrait with Equilibrium Points", labels = ["Rabbits", "Wolves"])

Now, we need to handle a modified version of the Lotka-Volterra equations. These modified equations incorporate logistic growth fot the prey population.

 

 

restart

alpha := 1.0; beta := .1; g := 1.5; delta := 0.75e-1; k := 100; ode1 := diff(x(t), t) = alpha*x(t)*(1-x(t)/k)-beta*x(t)*y(t); ode2 := diff(y(t), t) = delta*x(t)*y(t)-g*y(t); eq1 := alpha*x*(1-x/k)-beta*x*y = 0; eq2 := delta*x*y-g*y = 0; equilibria := solve({eq1, eq2}, {x, y}); print("Equilibrium Points: ", equilibria); initial_conditions := x(0) = 40, y(0) = 9; sol := dsolve({ode1, ode2, initial_conditions}, {x(t), y(t)}, numeric); eq_points := [seq([rhs(eq[1]), rhs(eq[2])], `in`(eq, equilibria))]

[[0., 0.], [100., 0.], [20., 8.]]

plots[odeplot](sol, [[t, x(t)], [t, y(t)]], t = 0 .. 100, legend = ["Rabbits", "Wolves"], title = "Prey-Predator Dynamics", labels = ["Time", "Population"])

NULL

plots:-odeplot(sol, [[x(t), y(t)]], 0 .. 50, color = "blue"); equilibrium_plot := plots:-pointplot(eq_points, color = "red", symbol = solidcircle, symbolsize = 15); plots:-display([plots:-odeplot(sol, [[x(t), y(t)]], 0 .. 50, color = "blue"), equilibrium_plot], title = "Phase Portrait with Equilibrium Points", labels = ["Rabbits", "Wolves"])

NULL

 

Download predator_prey2.mw

solving 2 equations...

dharr 8200
August 05 2024
1 2

In your first section of code, after the first solve, you assign x[2]^2/3 to x[1]. So now fx is (x[2]^2/3)^2 - 3*x[2] = 0, so a quartic equation in x[2]. When you solve it you get the 4 roots of this quartic. So this is correct.

A better approach is to just solve the two equations simultaneously, i.e., let Maple do the work for you.

sol := solve({fx, fy},{x[1],x[2]})


you could then use assign to set the values of x[1] and x[2], or better, use the solution you want with eval to work out other things involving x[1] and x[2].
(For future reference, please upload your worksheet contents and leave the link, so we can directly use it without retyping it in.)

> 

restart
> 

F4 := x[1]^3+x[2]^3-9*x[1]*x[2]+27

x[1]^3+x[2]^3-9*x[1]*x[2]+27

> 

fx := (1/3)*(diff(F4, x[1])); fy := (1/3)*(diff(F4, x[2]))

x[1]^2-3*x[2]

x[2]^2-3*x[1]

[] are optional below, but allows you to, e.g, count solutions

> 

sol := [solve({fx, fy}, {x[1], x[2]}, explicit)]; nops(sol)

[{x[1] = 0, x[2] = 0}, {x[1] = 3, x[2] = 3}, {x[1] = -3/2-((3/2)*I)*3^(1/2), x[2] = -3/2+((3/2)*I)*3^(1/2)}, {x[1] = -3/2+((3/2)*I)*3^(1/2), x[2] = -3/2-((3/2)*I)*3^(1/2)}]

4

Example: value of F4 for the 4th solution

> 

eval(F4, sol[4]); simplify(%)

(-3/2+((3/2)*I)*3^(1/2))^3+(-3/2-((3/2)*I)*3^(1/2))^3-9*(-3/2+((3/2)*I)*3^(1/2))*(-3/2-((3/2)*I)*3^(1/2))+27

0

> 

NULL

Download solve.mw

 

InertForm...

dharr 8200
August 05 2024
1 3

I assume you are OK with the standard mathematical way of writing the limit. Then the following works.

InertForm:-Display(`%>`(%limit(X[i] %^ (rho = 0), alpha = infinity), 0), inert = false)

 

one sqrt...

dharr 8200
August 02 2024
3 0

Maple calculates THE sqrt (sqrt is a function); sqrt(4) does not give +/-2. From the help page "sqrt(x) represents the "principal square root", defined by the formula sqrt(x) = exp(1/2 * ln(x))". 

coeffs or solve,identity...

dharr 8200
August 02 2024
0 0

coeffs finds all the coefficients, and then they can be used in solve (by default solve assumes they are equal to zero). But solve(identity(,...) automates this process. There is no solution for eq4, but for eq1 if you do not specify q[0], then there is a trivial solution.

Solve_for_coefficients2.mw

bug...

dharr 8200
August 01 2024
1 0

I'd say it is a bug. A workaround if to use evalf

sol:=[-1/4*x^2, (-1/2-1/2*(-3)^(1/2))^2+(-1/2-1/2*(-3)^(1/2))*x, (-1/2+1/2*(-3)^(1/2))^2+(-1/2+1/2*(-3)^(1/2))*x];
plot(evalf(sol),legend=sol)

 

First 14 15 16 17 18 19 20 Last Page 16 of 81