Advanced Search

dharr

Dr. David Harrington

4444 Reputation

20 Badges

18 years, 280 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 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 1990.

MaplePrimes Activity


These are answers submitted by dharr

Using dsolve with NLPSolve...

dharr 4444
May 31 2023
1 7

Here's how I do least squares fits with numerical solutions from dsolve. I bypass the parameters facility of dsolve. I've only tested this on the derivative-free nonlinearsimplex method [but see below]. I use the matrix method call to NLPSolve; for the procedure method it seemed I had to hard-code the number of parameters in the internal procedure SS.

Edit: it seems to work with sqp, but the initial point had to be closer.

> 

restart;

Test DE with parameters in DE and initial conditions

> 

de:={diff(y(t),t,t)=-a^2*y(t),D(y)(0)=b*a,y(0)=c};

{diff(diff(y(t), t), t) = -a^2*y(t), y(0) = c, (D(y))(0) = b*a}

Analytical answer

> 

ans:=rhs(dsolve(de,y(t)));

b*sin(a*t)+c*cos(a*t)

Generate signal for some test parameters

> 

params:={a=1.1,b=2.,c=3.};
sig:=unapply(eval(ans,params),t);
tt:=Vector([seq(t,t=0..6,0.1)]):
yy:=map(sig,tt):

{a = 1.1, b = 2., c = 3.}

proc (t) options operator, arrow; 2.*sin(1.1*t)+3.*cos(1.1*t) end proc

> 

#plot(tt,yy,style=point);

Use matrix method call to NLPSolve

> 

defit:=proc(tims::~Vector,expt::~Vector,des::set,inivarsarg::list,{tvar::symbol:='t',demethod:='rkf45',fitmethod:='nonlinearsimplex'})
  local ans,params,gparams,vars,inivals,SS;
  vars:=[indets(des,function(symbol))[]];
  gparams:=convert(indets(des,'name') minus {tvar},list);
  inivals:=eval(gparams,inivarsarg);
  SS:=proc(V)
    local ans;
    if not type(V,'Vector'(numeric)) then return 'procname'(args) end if;  
    ans:=dsolve(eval(des,gparams=~convert(V,list)),vars,numeric,output=Array(tims),method=demethod);
    add((expt-ans[2,1][..,2])^~2);
  end proc:
  Optimization:-NLPSolve(nops(inivarsarg),SS,initialpoint=Vector(inivals),method=fitmethod);
end proc:
> 

defit(tt,yy,de,[a=1.5,b=2.5,c=3.3],tvar=t,fitmethod=sqp,demethod=rkf45);

[1.33021633758686053*10^(-11), Vector(3, {(1) = 1.0999999857061296, (2) = 1.999999370723821, (3) = 2.9999993236548868})]

> 

``

Download defit.mw

assign...

dharr 4444
May 27 2023
0 0
> 

restart;

You already have the names and values in tables, I assume something like this:

> 

varnames[1]:=varname1;
varnames[2]:=varname2;
values[1]:=0.600;
values[2]:=123;

varname1

varname2

.600

123

Now assign them

> 

assign(seq(varnames[i]=values[i],i=1..numelems(varnames)));
> 

varname1;
varname2;

.600

123

> 

NULL

Download assign.mw

You can do the assign as you read a line by putting a line such as "splog=0.56" into a variable str, and then use

assign(parse(str));

convex hull...

dharr 4444
May 21 2023
5 1

In this case, the points of the polyhedron and its convex hull are the same, so ConvexHull in the ComputationalGeometry package works

> 

restart
> 

with(ComputationalGeometry)
> 

pts := [[1, 1, 0], [-1, 1, 0], [-1, -1, 0], [-1/2, 0, 1], [1/2, -1/2, 0]]

[[1, 1, 0], [-1, 1, 0], [-1, -1, 0], [-1/2, 0, 1], [1/2, -1/2, 0]]

> 

polyhedron := ConvexHull(pts)

[[1, 4, 2], [4, 3, 2], [5, 4, 1], [4, 5, 3], [3, 5, 2], [5, 1, 2]]

> 

plots:-display(map(proc (tri) options operator, arrow; plottools:-polygon(pts[tri], transparency = .7, color = red) end proc, polyhedron))

> 

ConvexHull(pts, output = volume)

1.

> 

NULL

Download volume.mw

x range, numpoints...

dharr 4444
May 20 2023
1 0

If you work out the exact x range for the spacecurve it works, e,g, for the circle it is just x=-3..3. For the ellipse see the attached. Alternatively, you can just use numpoints=10000

SectionsConiquesTest.mw.

different T...

dharr 4444
May 17 2023
2 8

The T you used in ConvertIn(T) is not the same as the internal one. The easiest way to get around this is always to provide an irreducible polyomial to GF, then the variable you use there can be used.

> 

restart;
> 

p:=Nextprime(T^4,T) mod 2;
G:=GF(2,4,p);
b:=G:-ConvertIn(T^2+T);
G:-ConvertOut(b);

T^4+T+1

_m1832374234528

modp1(ConvertIn(T^2+T, T), 2)

T^2+T

> 

NULL

Download GF.mw

CUDA...

dharr 4444
May 16 2023
2 29

CUDA only helps speed up numerical floating point matrix multilpication, which will not help your symbolic solution of polynomial equations. See ?CUDA,supported_routines

one way...

dharr 4444
May 13 2023
0 0
> 

restart
> 

a[1] := -d[l]*d[p]+p[2]*rho*A*r[1]*d[i]/d[u]-mu*d[i]-mu*d[p]-d[i]*d[p]-d[i]*d[l]

-d[l]*d[p]+p[2]*rho*A*r[1]*d[i]/d[u]-mu*d[i]-mu*d[p]-d[i]*d[p]-d[i]*d[l]

Make this an equation

> 

eq := R[0] = (r[1]*p[2]/d[p]+mu*r[1]*p[1]/((d[l]+mu)*d[p]))*A*rho/d[u]

R[0] = (r[1]*p[2]/d[p]+mu*r[1]*p[1]/((d[l]+mu)*d[p]))*A*rho/d[u]

Suppose we want to substitute, eliminating A (other choices are possible)

> 

A = solve(eq, A); subs(%, a[1]); simplify(%)

A = R[0]*(d[l]+mu)*d[p]*d[u]/(r[1]*(mu*p[1]+mu*p[2]+d[l]*p[2])*rho)

-d[l]*d[p]+p[2]*R[0]*(d[l]+mu)*d[p]*d[i]/(mu*p[1]+mu*p[2]+d[l]*p[2])-mu*d[i]-mu*d[p]-d[i]*d[p]-d[i]*d[l]

-d[l]*d[p]+p[2]*R[0]*(d[l]+mu)*d[p]*d[i]/((d[l]+mu)*p[2]+mu*p[1])-mu*d[i]-mu*d[p]-d[i]*d[p]-d[i]*d[l]

NULL

Download subs.mw

indices...

dharr 4444
May 13 2023
2 1

You need A[1] and A[2], not A__1 and A__2 (they look similar in the output). A__1 is just a name that appears as A with a subscript, but it is not related to the Array A.

plots.mw

units; linear and nonlinear regression...

dharr 4444
May 10 2023
3 2

To get the values right, you need R in kOhms, C in mF. The original formula is nonlinear and the sumX etc are for linear regression. If we take logs first them the linear regression works. Otherwise we can use Maple's built-in nonlinear fitting routine.

> 

restart

Use U in volts, R in kOhms, C in mF so RC is in seconds, U/R is mA

> 

U := 100

100

> 

R := 4.900

4.900

> 

C := 20*(1/1000)

1/50

> 

"v(t):=U/(R)*(e)^(-(t/(R*C)))"

proc (t) options operator, arrow, function_assign; U*exp(-t/(R*C))/R end proc

> 

v(.1)

7.356077314

``

> 

data := LinearAlgebra:-Transpose(`<,>`(`<|>`(.1, .2, .3, .4, .5), `<|>`(7.36, 2.7, .99, .37, .13)))

Matrix(%id = 36893490582268737284)

Plot is nonlinear

> 

plot(data, style = point)

Take logs to make it linear

> 

X := `<,>`(seq(.1 .. .5, .1)); Y := map(proc (t) options operator, arrow; ln(v(t)) end proc, X); lndata := `<|>`(X, Y)

Vector(5, {(1) = .1, (2) = .2, (3) = .3, (4) = .4, (5) = .5})

Vector(5, {(1) = 1.995526817, (2) = .9751186554, (3) = -0.4528950891e-1, (4) = -1.065697672, (5) = -2.086105834})

Matrix(%id = 36893490582268711868)

Now it is linear

> 

plot(lndata, style = point)

> 

sumX := add(X); sumY := add(Y); sumX2 := add(x^2, `in`(x, X)); sumXY := X.Y

1.5

-.226447543

.55

-1.08834242569299988

> 

eq1 := numelems(X)*a+b*sumX = sumY; eq2 := a*sumX+b*sumX2 = sumXY

5*a+1.5*b = -.226447543

1.5*a+.55*b = -1.08834242569299988

> 

sol := solve({eq1, eq2}, {a, b})

{a = 3.015934980, b = -10.20408163}

Slope of ln(v(t) = ln(U/R) - t/(R*C) is -1/RC, intercept is ln(U/R)

> 

-1/(R*C)

-10.20408164

> 

ln(U/R)

3.015934981

Or we could use nonlinear fitting to fit to the original equation

> 

C := 'C'; R := 'R'

C

R

> 

v(t)

100*exp(-t/(R*C))/R

> 

Statistics:-Fit(v(t), data, t, output = parametervalues)

[C = HFloat(0.02000589855860291), R = HFloat(4.986433379661445)]

NULL

Download regression.mw

package loading...

dharr 4444
May 09 2023
3 2

The only thing that needs to be done is to set libname correctly. The simplest way to get started is to open the examples.mws file from where it is in the unzipped directory. Then insert a line after restart which has

libname:=libname,interface(worksheetdir);

Then the package should load.

All that needs to happen is that libname points to wherever the maple.ind,maple.hdb and maple.lib are located (I would rename them FourierTrigSeries.lib etc to make their contents clearer). The above works because these files are in the same directory as the worksheet.

This libname command can later be put into a maple.ini file to make it automatic, as described in the readme file.

sqrt for normal matrices...

dharr 4444
May 07 2023
4 5

Maple uses an algorithm that works for all matrices, but the square root of normal matrices can be done more simply.
[Edit: improved version of UDiag]

Matrices in startup code.

> 

restart
> 

with(LinearAlgebra)

M__4 is symmetric (therefore normal), and therefore we can do a matrix function on its eigenvalues

> 

Equal(M__4, LinearAlgebra:-Transpose(M__4))

true

Unitary diagonalization of normal matrix, A = U.Lambda.U^*.  Normal matrices, i.e, A.A^*=A^*.A have no missing eigenvectors, i.e., geometric and algebraic multiplicities the same.

The Eigenvectors command produces eigenvectors that are already orthogonal if the eigenvalues are distinct, so they need to be normalized, and orthogonalized if they belong to the same eigenvalue.

> 

UDiag:=proc(A::Matrix(square))
  local evs,evals,evecs,U,i,j,k,g,n;
  uses LinearAlgebra;
  n:=RowDimension(A);
  evs:=Eigenvectors(A,output='list');
  evals:=Vector(n);evecs:=Matrix(n,n);
  j:=0;
  for i in evs do
    if i[2] = 1 then
      j:=j+1;
      evals[j]:=i[1];
      evecs[..,j]:=Normalize(i[3][],'Euclidean');
    else
      g:=GramSchmidt(i[3],conjugate=false,normalized=true);
      for k to i[2] do
        j:=j+1;
        evals[j]:=i[1];
        evecs[..,j]:=g[k];
      end do;
    end if;
  end do;
  DiagonalMatrix(evals),evecs;
end proc:
> 

Lambda, U := UDiag(M__4)
> 

sqrtM__4 := map(simplify, U.map(sqrt, Lambda).LinearAlgebra:-HermitianTranspose(U))
> 

Equal(map(simplify, sqrtM__4.sqrtM__4), M__4)

true

> 

Lambda, U := UDiag(M__1); sqrtM__1 := map(simplify, U.map(sqrt, Lambda).LinearAlgebra:-HermitianTranspose(U)); Equal(map(simplify, sqrtM__1.sqrtM__1), M__1)

true

> 

Lambda, U := UDiag(M__2); sqrtM__2 := map(simplify, U.map(sqrt, Lambda).LinearAlgebra:-HermitianTranspose(U)); Equal(map(simplify, sqrtM__2.sqrtM__2), M__2)

true

> 

Lambda, U := UDiag(M__3); sqrtM__3 := map(simplify, U.map(sqrt, Lambda).LinearAlgebra:-HermitianTranspose(U)); Equal(map(simplify, sqrtM__3.sqrtM__3), M__3)

true

> 

NULL

Download mpower.mw

tables...

dharr 4444
May 03 2023
0 1

Indexed notation like cov[x]:=newval; would be the usual way to modify an entry. I don't see anything wrong with your code; I'd use eval(cov) over op(cov) for the return value for a slight improvement in readability. To reset a table, just set it to a new blank table: cov:=table();

In Maple 2023, cov := table([seq(x = x, x = X)]) could be cov := table(X =~ X)

and cov := table([seq(x = NULL, x = X)]) could be cov:=table(X =~ NULL)

RootOf, fsolve stuff...

dharr 4444
April 30 2023
1 1

Negative values (as opposed to index=negative value) are just estimates of roots. index and allvalues work best for polynomials. I do it also by fsolve, by you can also use NextZero for systematically finding roots.

I don't know another way to add "index=number", other than convert(r1,RootOf,form=index) which gives index=1 (for a polynomial).

> 

restart
> 

arg := _Z*cos(_Z)-sqrt(sin(_Z)^2)

_Z*cos(_Z)-(sin(_Z)^2)^(1/2)

> 

r0 := RootOf(arg)

RootOf(_Z*cos(_Z)-(sin(_Z)^2)^(1/2))

> 

inter := -6 .. 6

-6 .. 6

> 

plot(arg, _Z = inter, size = [300, 300])

Get any one root

> 

evalf(r0); fsolve(arg); fsolve(arg, _Z = inter)

0.

0.

-2.028757838

If you are willing to look at the plot and do hand estimates then RootOf (or fsolve) works. Save this input line first, and then evaluate.

> 

vals := 'evalf(RootOf(arg, -4)), evalf(RootOf(arg, -2)), evalf(RootOf(arg, 0)), evalf(RootOf(arg, 5))'

evalf(RootOf(arg, -4)), evalf(RootOf(arg, -2)), evalf(RootOf(arg, 0)), evalf(RootOf(arg, 5))

> 

vals

-4.493409458, -2.028757838, 0., 4.913180439

Or intervals are probably more reliable.

> 

evalf(RootOf(arg, -5 .. -4)), evalf(RootOf(arg, -3 .. -1)), evalf(RootOf(arg, -1 .. 1)), evalf(RootOf(arg, 4 .. 5))

-4.493409458, -2.028757838, 0., 4.913180439

Conversions are possible - the strange result here is because index is really  only for polynomials

> 

r1 := RootOf(arg, -4); convert(r1, RootOf, form = index); evalf(%)

RootOf(_Z*cos(_Z)-(sin(_Z)^2)^(1/2), -4)

RootOf(_Z*cos(_Z)-(eval(RootOf(_Z^2-E, index = 1), {E = sin(_Z)^2})), -4)

-4.493409458

Not sure whats happening in the next two

> 

convert(r1, 'RootOf', form = 'interval')

Error, (in type/algext) too many levels of recursion

> 

convert(RootOf(arg, -5 .. -4), RootOf, form = numeric)

Error, (in type/algext) too many levels of recursion

Using fsolve and avoid

> 

a1 := fsolve(arg, _Z = inter); a2 := fsolve(arg, _Z, avoid = {_Z = a1}); a3 := fsolve(arg, _Z, avoid = {_Z = a1, _Z = a2}); a4 := fsolve(arg, _Z, avoid = {_Z = a1, _Z = a2, _Z = a3})

-2.028757838

0.

-4.493409458

4.913180439

> 

NULL

Download RootOf.mw

extra spaces...

dharr 4444
April 30 2023
1 8

In future it would be helpful if you uploaded your worksheet with the green up-arrow in the Mapleprimes editor. One problem is that your code shows diff*(S(t), t) not diff(S(t), t) because you have a space between diff and (, which is interpreted as multiplication. (also in other places).

triangulations...

dharr 4444
April 29 2023
4 0

Here is a solution, based on iterating through the related binary trees - see https://en.wikipedia.org/wiki/Catalan_number

> 

restart

Consider a regular polygon with m = n+2 vertices. We want all non-crossing triangulations of the vertices, see https://en.wikipedia.org/wiki/Catalan_number 
These are in 1:1 correspondence with the binary trees with n+1 leaves or n internal nodes, or with n deep nested parentheses.
Consider the example of a pentagon.

> 

with(GraphTheory)
> 

m := 5; n := m-2

5

3

The edges of the pentagon are labelled A,B,C,D,1. 1 corresponds to the first level vertex of the binary tree, and the others to the leaves. Edges 2..n will be internal edges added, which correspond to the internal vertices of the binary tree.

> 

leaves := [seq({i, i+1}, i = 1 .. m-1)]; edges := {leaves[], {1, m}}; G := Graph(edges)

[{1, 2}, {2, 3}, {3, 4}, {4, 5}]

{{1, 2}, {1, 5}, {2, 3}, {3, 4}, {4, 5}}

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {2, 5}, (2) = {1, 3}, (3) = {2, 4}, (4) = {3, 5}, (5) = {1, 4}}), `GRAPHLN/table/1`, 0)

The P representation (see ?Iterator:-Trees) for binary trees/nested parentheses gives the internal tree vertices, which are interspersed with the leafs to give a nested parentheses representation,

e.g. [1,3,2] means A 1 B 3 C 2 D or  A1((B3C)2D) or just A((BC)D), where the 3 indicates the first pairing and 2 the second and 1 the last.
Starting with 3, we generate internal edge 3 as the third side of the triangle B3C, then internal edge 2 is the third side of the triangle 32D, then the last triangle is 21A

> 

p := Array([1, 3, 2])

Array(%id = 36893489648404452884)

The following (inelegant and inefficient) code triangulates this case, contracting the parentheses in turn

> 

plist := convert(p, list);
edgelist := leaves;
for i from n by -1 to 2 do
    member(i, plist, 'j');
    print(i, j);
    newedge := symmdiff(edgelist[j], edgelist[j + 1]);
    AddEdge(G, newedge);
    edgelist := [edgelist[1 .. j - 1][], newedge, edgelist[j + 2 .. -1][]];
    plist := subsop(j = NULL, plist);
end do;

[1, 3, 2]

[{1, 2}, {2, 3}, {3, 4}, {4, 5}]

true

3, 2

{2, 4}

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {2, 5}, (2) = {1, 3, 4}, (3) = {2, 4}, (4) = {2, 3, 5}, (5) = {1, 4}}), `GRAPHLN/table/1`, 0)

[{1, 2}, {2, 4}, {4, 5}]

[1, 2]

true

2, 2

{2, 5}

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {2, 5}, (2) = {1, 3, 4, 5}, (3) = {2, 4}, (4) = {2, 3, 5}, (5) = {1, 2, 4}}), `GRAPHLN/table/1`, 0)

[{1, 2}, {2, 5}]

[1]

> 

DrawGraph(G, size = [200, 200])

So here is a procedure that iterates though all cases

> 

triangulations:=proc(m::posint)
  local n,leaves,edges,i,t,j,G,T,p,plist,newedge,graphs,edgelist,ng;
  uses GraphTheory;
  n:=m-2;
  leaves := [seq({i, i + 1}, i = 1 .. m - 1)];
  edges := {leaves[], {1, m}};
  T:=Iterator:-BinaryTrees(n);
  graphs:=table();
  ng:=0;
  for t in T do
    G := Graph(convert(edges, set));
    p := Iterator:-Trees:-Convert(t, 'LR', 'P');
    plist := convert(p, list);
    edgelist := leaves;
    for i from n by -1 to 2 do
      member(i, plist, 'j');
      newedge := symmdiff(edgelist[j], edgelist[j + 1]);
      AddEdge(G, newedge);
      edgelist := [edgelist[1 .. j - 1][], newedge, edgelist[j + 2 .. -1][]];
      plist := subsop(j = NULL, plist);
    end do;
    ng:=ng+1;
    graphs[ng]:=copy(G);
  end do;
  convert(graphs,list);
end proc:

And here are the 14 cases for m=6

> 

DrawGraph(triangulations(6))

 

 
 

 

 

> 

NULL

Download Catalan.mw

1 2 3 4 5 6 7 Last Page 1 of 51