dharr

Dr. David Harrington

6892 Reputation

21 Badges

20 years, 108 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 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 replies submitted by dharr

I assume grade means degree here? Some further details about typical coeffs, number of terms, degree etc would be helpful. For integer coeffs, dense or sparse number of terms, degree 10^5 seems feasible, but with both degrees 10^6, it took longer than I was willing to wait. (This was Maple 2024 on Windows 11.)
 

restart;

r:=randpoly(x, degree = 15);
p:=expand(r*randpoly(x, dense, degree=10^5)):
q:=expand(r*randpoly(x, dense, degree=10^5+20)):

87*x^10-56*x^9-62*x^3+97*x-73

CodeTools:-Usage(gcd(p,q));

memory used=70.30MiB, alloc change=35.88MiB, cpu time=4.72s, real time=19.49s, gc time=0ns

87*x^10-56*x^9-62*x^3+97*x-73

p:=expand(r*randpoly(x, degree=10^6)):
q:=expand(r*randpoly(x, degree=10^5)):

CodeTools:-Usage(gcd(p,q));

memory used=8.03MiB, alloc change=-29.05MiB, cpu time=1.06s, real time=3.67s, gc time=0ns

(87*x^10-56*x^9-62*x^3+97*x-73)*x^51516

 

Download gcd.mw

 

@MaPal93 I'm not sure how they are chosen. They seem to be points with integer coordinates on a grid, but I suppose if a region was too small for that then some rationals would be chosen. Note that for regions 2 and 3, the horizontal axis has no particular significance. I think the basic idea is to have any point for which evaluation proves that the equations and inequalities are satisfied; in that sense it doesn't matter where in the region it is.

@C_R With floating point limits, it is supposed to switch to numerical integration. From ?int

"If any of the integration limits of a definite integral are floating-point numbers (e.g. 0.0, 1e5 or an expression that evaluates to a float, such as exp(-0.1)), then int computes the integral using numerical methods if possible (see evalf/Int). Symbolic integration will be used if the limits are not floating-point numbers unless the numeric=true option is given."

But the alternate form also gave a numeric answer with 0..1,0..1 which is not expected, unless the floats in the integrand triggered this.

@mmcdara Very nice. Vote up. 

For axes=normal, the default, I don't know how to change this, but if you are OK with axes=boxed, then the labels are on the other side.

@C_R So even for polynomials allvalues will use radicals for cubics unless told not to. So to avoid radicals, and possible missed solutions, it should be possible to tell dsolve to use solve to get a RootOf. It turns out it is not so simple, but it can be done without as much done "by hand" as I did above.

restart

dsolve without method

ode:=diff(y(x), x) = (3*x - y(x) + 1)/(3*y(x) - x + 5);
ic:=y(0)=0;
dsolve({ode,ic});

diff(y(x), x) = (3*x-y(x)+1)/(3*y(x)-x+5)

y(0) = 0

y(x) = -(-(1/36)*(x+1)^2*((-324+12*(96*x^3+288*x^2+288*x+825)^(1/2))^(2/3)-24*x-24)^2/(-324+12*(96*x^3+288*x^2+288*x+825)^(1/2))^(2/3)-x^3-x^2+x+1)/(x+1)^2

According to the help, (?dsolve,details), dsolve is sensitive to _EnvExplicit. so I expected this to work

_EnvExplicit:=false;

false

but no

dsolve({ode,ic});

y(x) = -(-(1/36)*(x+1)^2*((-324+12*(96*x^3+288*x^2+288*x+825)^(1/2))^(2/3)-24*x-24)^2/(-324+12*(96*x^3+288*x^2+288*x+825)^(1/2))^(2/3)-x^3-x^2+x+1)/(x+1)^2

But if we dsolve implicitly, now solve will use RootOf. It seems to find each solution 3 times, but it doesn't miss any.

dsolve({ode,ic},implicit);
solve(%,y(x));
plot([allvalues(%,implicit)],x=-10..10);

-(2/3)*ln(-(3+y(x)+x)/(x+1))-(1/3)*ln((-1-y(x)+x)/(x+1))-ln(x+1)+(2/3)*ln(3)+I*Pi = 0

-RootOf(9+(x^3+3*x^2+3*x+1)*_Z^9+(2*x^3+6*x^2+6*x+2)*_Z^6)^3*x-RootOf(9+(x^3+3*x^2+3*x+1)*_Z^9+(2*x^3+6*x^2+6*x+2)*_Z^6)^3-x-3

NULL

 

NULL

Download RootOf_stuff_2.mw

@Axel Vogt Yes, I noticed that early on, but it was unclear about what exactly the poster wanted. If the poster can explain the practical range of x that they want the integral for I might be willing to look more closely. I've had success with method = _do1akc in the past with oscillating integrals. But non-dimensionalization would be a good first step; I doubt reliable numerics can be found with the parameters ranging from 10^23 to 10^(-37). m is set to 1 but I'm not sure if this is just a test parameter. If it is later to have very different values that also would make a big difference since it multiplies x.

We are still missing a value for x.

The first integral is problematic, since the integrand

(-epsilon^2/(4*a^2) - tau*epsilon^3)*(cos(epsilon*x) - 1)/epsilon^2

goes to infinity as epsilon goes to infinity because of the tau* epsilon^3. Is the first factor perhaps intended to be within an exp?

@Carl Love OK, thanks for that.

If you reorganize the dataframe so the your Data are the row labels (which seems logical to me but may not be what you want), then you can just select the rows and columns by indexing:

TestData[["LC3", "LC4", "LC5", "LC6", "LC7", "LC8", "LC9"], Load];

dataframe.mw

@C_R I think Maple intends polar to serve that purpose - you can do arithmetic with it, but you need evalc() or simplify() to get back to the Re/Im form.

restart

a := convert(-2, polar)

polar(2, Pi)

b := convert(7, polar)

polar(7, 0)

a*b; simplify(%)

polar(2, Pi)*polar(7, 0)

-14

a+b; simplify(%)

polar(2, Pi)+polar(7, 0)

5

c := convert(2+I, polar)

polar(5^(1/2), arctan(1/2))

a+c; simplify(%)

polar(2, Pi)+polar(5^(1/2), arctan(1/2))

I

d := convert(3+I, polar)

polar(10^(1/2), arctan(1/3))

a+d; simplify(%); evalc(%)

polar(2, Pi)+polar(10^(1/2), arctan(1/3))

-2+polar(10^(1/2), arctan(1/3))

1+I

NULL

Download polar.mw

@Carl Love Generally true of course, but the statement needs some qualification.

eq1 := x*y = 1

x*y = 1

ans := solve(eq1, {x, y})

{x = 1/y, y = y}

Sample solution

map(proc (q) options operator, arrow; lhs(q) = eval(rhs(q), y = 2) end proc, ans)

{x = 1/2, y = 2}

NULL

Download solve.mw

I am unable to reproduce this (Maple 2024, Windows 11). After the reformatting thw worksheet runs correctly.

cost_comparison_-_liquid_(v01MP).mw

@dharr Of course you can just write X = A.L, with A = X.L^(-1). Still not sure what the rules are.

@WA573 

restart

with(LinearAlgebra)

L := Matrix([[(1+I*w*sqrt((1-cos(k))/(1+cos(k)))/lambda)*exp(-(1/2)*k), I*rho*(exp((1/2)*k)-exp(-(1/2)*k))/lambda], [I*rho*(exp((-k)*(1/2))-exp((1/2)*k))/lambda, (1-I*w*sqrt((1-cos(k))/(1+cos(k)))/lambda)*exp(-(1/2)*k)]]); tL := simplify(Trace(L))

Matrix(%id = 36893490171766514140)

2*exp(-(1/2)*k)

X := Matrix([[-I*lambda*(1/2)-(1/2)*w, -rho], [rho, I*lambda*(1/2)+(1/2)*w]]); Trace(X)

Matrix(%id = 36893490171766521972)

0

X has a pair of eigenvalues = +/-mu; their average is at the origin. To shift the eigenvalues of L so their average is at the origin, we add a matrix S

S := -(1/2)*tL*IdentityMatrix(2); Q := S+L; simplify(Trace(Q))

Matrix(2, 2, {(1, 1) = -exp(-(1/2)*k), (1, 2) = 0, (2, 1) = 0, (2, 2) = -exp(-(1/2)*k)})

Matrix(%id = 36893490171884120892)

0

Now a rotation/contraction in the complex plane by multiplying by the ratio of the eigenvalues b will make bQ and X similar

b := Eigenvalues(X)[1]/Eigenvalues(Q)[1]; bQ := b*Q

(1/2)*((2*I)*lambda*w-lambda^2-4*rho^2+w^2)^(1/2)*(1+cos(k))*exp(-(1/2)*k)*lambda/((1+cos(k))*((exp(-(1/2)*k))^4*cos(k)*rho^2+(exp(-(1/2)*k))^4*w^2*cos(k)+(exp(-(1/2)*k))^4*rho^2-(exp(-(1/2)*k))^4*w^2-2*(exp(-(1/2)*k))^2*cos(k)*rho^2-2*(exp(-(1/2)*k))^2*rho^2+cos(k)*rho^2+rho^2))^(1/2)

Eigenvalues(bQ, output = list); Eigenvalues(X, output = list)

[(1/2)*((2*I)*lambda*w-lambda^2-4*rho^2+w^2)^(1/2), -(1/2)*((2*I)*lambda*w-lambda^2-4*rho^2+w^2)^(1/2)]

[(1/2)*((2*I)*lambda*w-lambda^2-4*rho^2+w^2)^(1/2), -(1/2)*((2*I)*lambda*w-lambda^2-4*rho^2+w^2)^(1/2)]

Now bQ and X are similar. So what Matrix relates them?

evalsbQ, evecsbQ := Eigenvectors(bQ)

evalsX, evecsX := Eigenvectors(X)

So we have 1/evecsX.X.evecsX = 1/evecsbQ.bQ.evecsbQ or X = evecsX.(1/evexcbQ).bQ.evecsbQ.(1/evecsX)

B := evecsX.(1/evecsbQ); simplify(X-B__.bQ.(1/B))

Matrix(%id = 36893490171955216380)

So putting it all together we have X = b*B.Q.(1/B) and b*B.Q.(1/B) = b*B.(S+L).(1/B) and b*B.(S+L).(1/B) = -b.exp(-(1/2)*k).Identity+b*B.L.(1/B)

X__L := -b*exp(-(1/2)*k)*IdentityMatrix(2)+b*B.L.(1/B)

simplify(X-X__L)

Matrix(%id = 36893490171909539588)

NULL

NULL

NULL

NULL

 

NULL

Download XintoL4.mw

Earlier less streamlined version:

Download XintoL3.mw

First 8 9 10 11 12 13 14 Last Page 10 of 71