## dharr

Dr. David Harrington

## 4444 Reputation

18 years, 280 days
University of Victoria
Professor or university staff

## Social Networks and Content at Maplesoft.com

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.

## Using dsolve with NLPSolve...

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};

 > ans:=rhs(dsolve(de,y(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):

 > #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);

 >

## assign...

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

Now assign them

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

 >

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

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

 >
 >
 >

 >

 >

 >

 >

## x range, numpoints...

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

## different T...

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);

 >

## CUDA...

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

 >
 >

Make this an equation

 >

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

 >

## indices...

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

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.

 >

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

 >

 >

 >

 >

 >

 >

Plot is nonlinear

 >

Take logs to make it linear

 >

Now it is linear

 >

 >

 >

 >

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

 >

 >

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

 >

 >

 >

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);`

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

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.

 >
 >

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

 >

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

 >

 >

 >

 >

## tables...

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

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).

 >
 >

 >

 >

 >

Get any one root

 >

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.

 >

 >

Or intervals are probably more reliable.

 >

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

 >

Not sure whats happening in the next two

 >

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

 >

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

Using fsolve and avoid

 >

 >

## extra spaces...

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

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

 >

Consider a regular polygon with  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  leaves or  internal nodes, or with  deep nested parentheses.
Consider the example of a pentagon.

 >
 >

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.

 >

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

 >

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;

 >

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

 >

 >