 Obtain the tri-stimulus XYZ values from the CIE Color matching functions.

 Show the gamut of maximum chroma for the standard observer model with a D65 Illuminant.

 Approximate the white point of a Planckian source and compare to D65.

 Translate the maximum chroma gamut in xy to Lab (CIE L*a*b*) for perceived gamut (Violet and Magenta come together)

 Map the RGB color cube of fully saturated color into Lab and compare to perceivable colors.

10/6/15  Initial Document

•12/28/15 Improve RGB gamut with more data points: Procedures added for RGB to Lab: Wavlength Colors now based on CIEDE2000 model for Lab.                   

 Here is the latest version of this document, the MSL_data must be in a directory set in the mw file;

MSL_data.xlsx    Vision_RGB_Gamut.mw

In this paper we will demonstrate the many differences of implementation in the modeling of mechanical systems using embedded components through Maplesoft. The mechanical systems are used for different tasks and therefore have different structure in its design; as to the nature of the used functional elements placed on them, they vary greatly. This diversity is reflected in approaches and practices in modeling.

The following cases focus on mechanical components of the units manufacturing and processing machines. We can generate graphs for analysis using different dynamic pair ametros; all in real-time considerations in its manufacturing costs from the equations of conservation of energy.
Therefore modeling with Maplesoft ensures the smooth optimum performance in mechanical systems, highlighting the sustainability criteria for other areas of engineering.

XXXIII_Coloquio_SMP_2015.pdf

XXXIII_Coloquio_UNASAM_2015.mw

(in spanish)

L.AraujoC.

Apparently inconsistent behaviour of the BesselJ() function.

Examples: BesselJ(-3, 0)  ... gives 0 (correct)

but BesselJ(-3.0, 0), BesselJ(-3, 0.0)  and BesselJ(-3, 0.0) all give Float(infinity) (wrong! - should be 0.0)

The problem seems to occur for all negative integer values of the first argument (the order) when the second argument is 0 or 0.0.

Download Physics[Factor].mw

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

Maple's dsolve numeric can solve delay ODEs and DAEs as of Maple 18. However, if I am not wrong, it cannot solve delay equations with a time dependent history. In this post I show two examples.

Example 1:

y1(t) and y2(t) with time dependent history. Use of piecewise helps this problem to be solved efficiently. Hopefully Maple will add history soon in its capability.

Example 2: 

This is a very a complicated stiff problem from immunology. As of now, I believe only Maple can solve this (other than RADAR5 from Prof. Hairer). Details and plots are posted in the attached code.

Let me know if any one has a delay problem that needs to be solved. I have tested many delay problems in Maple (they work fine). The attached examples required addtional tweaking, hence the post.

I want to take this opportunity to congratulate and thank Maple's dsolve numeric/delay solvers for their fantastic job. Maple is world leader not because of example1, but because of its ability to solve example 2.

Download delayimmunetopost.mws

LL_104)_NASDAQ.mw
Portfolio_Optimization.txt

Portfolio Optimization with Google Spreadsheet and Maple
 

I will in this post show how to manage data and do portfolio optimization in Maple by using google spreadsheet.

You can either use a direct link to the data:

https://docs.google.com/spreadsheets/d/1L5-yUB0EWeBdJNMdELKBRmBQ1JJ0QymrtDLkVhHCVn8/pub?gid=649021574&single=true&output=csv

or you can set up your own google spreadsheet. If you choice to set up your own spreedsheet follow the below road map:

1) select which market you want to follow:

NASDAQ

http://www.nasdaq.com/screening/companies-by-industry.aspx?exchange=NASDAQ&render=download

NYSE

http://www.nasdaq.com/screening/companies-by-industry.aspx?exchange=NYSE&render=download

AMEX

http://www.nasdaq.com/screening/companies-by-industry.aspx?exchange=AMEX&render=download

2) Create a new google spreadsheet and name two sheets Blad1 and Panel. In the first cell of Blad1 you put the formula:

=IMPORTDATA("http://www.nasdaq.com/screening/companies-by-industry.aspx?exchange=NASDAQ&render=download")

you need to change the url to match your selection in 1).

3) In the first cell of Panel you put the name "Ticker" and then you copy all the ticker names from Blad1.

4) In the script editor you put in the below java script code:

function PanelCreation_Stock() 

{
var ss = SpreadsheetApp.getActiveSpreadsheet();
var sourceSheet = ss.getSheetByName("Blad1");
var dstSheet = ss.getSheetByName("Panel");
var curDat = new Date();
var day1 = curDat.getDay();
if(day1 == 0 || day1 == 1)
{
return;
}
var lCol = dstSheet.getLastColumn();
var srcdate = dstSheet.getRange(1, 1, 1, lCol).getValues();

for(var k=1;k<=srcdate[0].length-1;k++)
{
if(Utilities.formatDate(srcdate[0][k],"GMT", "dd-MMM-yy") == Utilities.formatDate(curDat,"GMT", "dd-MMM-yy"))
{
return;
}
}
var snRows = sourceSheet.getLastRow();
var dnRows = dstSheet.getLastRow();

var srcStock = sourceSheet.getRange("A2:A" + snRows).getValues();
var srcLastSale = sourceSheet.getRange("C2:C" + snRows).getValues();

var dstStock = dstSheet.getRange("A2:A" + dnRows).getValues();
var dstLastSale = dstSheet.getRange("Z2:Z" + dnRows).getValues();

for(var j=0;j<dnRows-1;j++)
{
dstLastSale[j][0]="n/a";
}
var flag = "true";
var foundStock;
for(var i=0;i<snRows-1;i++) //snRows
{
var sStockVal = srcStock[i][0];

//var foundStock = ArrayLib.indexOf(dstStock,0, sStockVal);

flag="false";
for(var j=0;j<dnRows-1;j++)
{
if(dstStock[j][0].toString().toUpperCase() == srcStock[i][0].toString().toUpperCase())
{
flag = "true";
foundStock = j;
break;
}
}
if(flag=="true")
{
dstLastSale[foundStock][0] = srcLastSale[i][0];
}
else
{
var dnRows1 = dstSheet.getLastRow()+1;
dstSheet.getRange("A" + dnRows1).setValue(srcStock[i][0]);
dstSheet.getRange(dnRows1,lCol+1,1,1).setValue(srcLastSale[i][0]);
for(var k=2;k<=lCol;k++)
{
if(dstSheet.getRange(dnRows1, k).getValue()=="")
{
dstSheet.getRange(dnRows1, k).setValue("n/a");
}
}
}
}
dstSheet.getRange(1,lCol+1).setValue(curDat);
dstSheet.getRange(2, lCol+1, dstLastSale.length, 1).setValues(dstLastSale);
}

 
5) Set it to run each day at 12:00. The code will save the new last sale price for monday to friday with one days lag.

Now we can move on to Maple.

In Maple run the following code to load the data:

X := proc (Url) local theDLL, URLDownloadToFile, myDirectory, myFile, Destination, DL;

theDLL := "C:\\WINDOWS\\SYSTEM32\\urlmon.dll";

URLDownloadToFile := define_external('URLDownloadToFileA', pCaller::(integer[4]), szURL::string, szFileName::string, dwReserved::(integer[4]), lpfnCB::(integer[4]), 'RETURN'::(integer[4]), LIB = theDLL);

if FileTools[Exists]("C:\\mydir") = true then FileTools:-RemoveDirectory("C:\\mydir", recurse = true, forceremove = true) else end if;

FileTools:-MakeDirectory("C:\\mydir");
myDirectory := "C:\\mydir";
myFile := "data1.csv";
Destination := cat(myDirectory, "\\", myFile);

DL := proc () local M;

URLDownloadToFile(0, Url, Destination, 0, 0);
M := ImportMatrix("C:\\mydir\\data1.csv", delimiter = ",", datatype = string);
M := Matrix(M, datatype = anything)

end proc;

return DL()

end proc:

data := X("https://docs.google.com/spreadsheets/d/1L5-yUB0EWeBdJNMdELKBRmBQ1JJ0QymrtDLkVhHCVn8/pub?gid=649021574&single=true&output=csv");
L := LinearAlgebra:-Transpose(data);

If you use your own spreadsheet you need to change the url to match that spreadsheet.
Select File -> Publish to the web in google spreadsheet

We can now run the portfolio optimization in Maple:

with(Statistics):
with(ListTools):
with(LinearAlgebra):
with(Optimization):
with(plots):

Nr, Nc := ArrayTools:-Size(L):
symb := L[1 .. 1, 2 .. Nc]:
LL := L[2 .. Nr, 2 .. Nc]:
Nr, Nc := ArrayTools:-Size(LL):

# Removing stocks with missing observations
for i to Nc do if Occurrences("n/a", convert(Column(LL, i), list)) >= 1 then AA[i] := i else AA[i] := 0 end if
end do;

DD := RemoveInRange([seq(AA[i], i = 1 .. Nc)], 0 .. 1):
symbb := DeleteColumn(symb, DD):
LLL := map(parse, DeleteColumn(LL, DD)):
Nr, Nc := ArrayTools:-Size(LLL):

# Calculate Return
for j to Nc do
for i from 2 to Nr do

r[i, j] := (LLL[i, j]-LLL[i-1, j])/LLL[i-1, j]

end do
end do;

RR := Matrix([seq([seq(r[i, j], j = 1 .. Nc)], i = 2 .. Nr)], datatype = float[8]);
n, nstock := ArrayTools:-Size(RR):

# Portfolio Optimization
W := Vector(nstock, symbol = w):
y := Vector(n, fill = 2, datatype = float[8]):
s1 := Optimization[LSSolve]([y, RR])[2];
Nr, Nc := ArrayTools:-Size(s1):

j := 0:
for i to Nr do if s1[i] <> 0 then j := j+1; ss1[j] := symbb[1, i] = s1[i] end if end do;

Vector(j, proc (i) options operator, arrow; ss1[i] end proc);
LineChart(s1);

Download PDEs_and_Boundary_Conditions.mw

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

Download SymbolicOrderDifferentiation.mw

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

A knight's tour is a sequence of moves of a knight on a chessboard such that the knight visits every square only once. This problem mentioned in  page of tasks still without  Maple implementation. 

The post presents the implementation of this task in Maple. Required parameter of the procedure (named  KnightTour)  is  address  - the address of the initial square in the algebraic notation. The second parameter  opt  is optional parameter:

1) if  opt is sequence  (by default) then  the procedure returns the sequence of moves of the knight in the usual algebraic notation,

2)  if  opt is diagram   then  the procedure returns the plot of moves of the knight and  sequentially numbers all the visited squares,

3) if  opt is animation  then  the procedure returns an animation of moves of the knight.

In the procedure is used a solution with maximum symmetry by George Jelliss, http://www.mayhematics.com/t/8f.htm

Code of the procedure:

KnightTour := proc(address::symbol, opt::symbol := sequence)

local L, n, L1, k, i, j, squares, border, chessboard, letters, digits, L2, L3, Tour, T, F;

uses ListTools, plottools, plots;

L := [a1, b3, d2, c4, a5, b7, d8, e6, d4, b5, c7, a8, b6, c8, a7, c6, b8, a6, b4, d5, e3, d1, b2, a4, c5, d7, f8, h7, f6, g8, h6, f7, h8, g6, e7, f5, h4, g2, e1, d3, e5, g4, f2, h1, g3, f1, h2, f3, g1, h3, g5, e4, d6, e8, g7, h5, f4, e2, c1, a2, c3, b1, a3, c2];

n := Search(address, L);

L1 := [L[n .. 64][], L[1 .. n-1][]];

if opt = sequence then return L1[] fi;

k := 0;

for i to 8 do

for j from `if`(type(i, odd), 1, 2) by 2 to 8 do

k := k+1;

squares[k] := polygon([[i-1/2, j-1/2], [i-1/2, j+1/2], [i+1/2, j+1/2], [i+1/2, j-1/2]], color = grey);

od;  od;

squares := convert(squares, list);

border := curve([[1/2, 1/2], [1/2, 17/2], [17/2, 17/2], [17/2, 1/2], [1/2, 1/2]], color = black, thickness = 4);

chessboard := display(squares, border);

letters := [a, b, c, d, e, f, g, h];

digits := [$ 1 .. 8];

L2 := convert~(L1, string);

L3 := subs(letters=~digits, map(t->[parse(t[1]), parse(t[2])], L2));

Tour := curve(L3, color = red, thickness = 3);

T := textplot([seq([op(L3[i]), i], i = 1 .. 64)], align = above, font = [times, bold, 16]);

if opt = diagram then return display(chessboard, Tour, T, axes = none) fi;

F := seq(display(chessboard, curve(L3[1 .. s], color = red, thickness = 3), textplot([seq([op(L3[i]), i], i = 1 .. s)], align = above, font = [times, bold, 16])), s = 1 .. 64);

display(seq(F[i]$5, i = 1 .. 64), insequence = true, axes = none);

end proc:

 Examples of use:

KnightTour(f3);

KnightTour(f3, diagram);

KnightTour(f3, animation);

 KnightTour.mw

Dear friends,

some time ago I shared a story here on the use of Maple to compute the cycle index of the induced action on the edges of an ordinary graph of the symmetric group permuting the vertices and the use of the Polya Enumeration Theorem to count non-isomorphic graphs by the number of edges. It can be found at the following Mapleprimes link.

I am writing today to alert you to another simple Maple program that is closely related and demonstrates Maple's capability to implement concepts from group theory and Polya enumeration. This link at Math.Stackexchange.com shows how to use the cycle index of the induced action by the symmetric group permuting vertices on the edges of a multigraph that includes loops to count set partitions of multisets containing two instances of n distinct types of items. The sequence that corresponds to these set partitions is OEIS A020555 where it is pointed out that we can equivalently count multigraphs with n labeled i.e. distinct edges where the vertices of the graph represent the multisets of the multiset partition and are connected by an edge k if the two instances of the value k are included in the sets represented by the two vertices that constitute the edge. The problem then reduces to a simple substitution into the aforementioned cycle index of a polynomial representing the set of labels on an edge including no labels on an edge that is not included.

This computation presents a remarkable simplicity while also implementing a non-trivial application of Polya counting. It is hoped that MaplePrimes users will enjoy reading this program, possibly profit from some of the techniques employed and be motivated to use Maple in their work on combinatorics problems.

Best regards,

Marko Riedel

Download MathematicalFunctionsSequences.mw

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

       Calculation of RSCR mechanism as a  solution to underdetermined system of nonlinear equations.  

RSCR.mw 

https://vk.com/doc242471809_376439263
https://vk.com/doc242471809_408704758

RCCC mechanism
https://vk.com/doc242471809_375452868

Solving DAEs in Maple

As I had mentioned in many posts, Maple cannot solve nonlinear DAEs. As of today (Maple 2015), given a system of index 1 DAE dy/dt = f(y,z); 0 = g(y,z), Maple “extends” g(y,z) to get dz/dt = g1(y,z). So, any given index 1 DAE is converted to a system of ODEs dy/dt = f, dz/dt = g1 with the constraint g = 0, before it solves. This is true for all the solvers in Maple despite the wrong claims in the help files. It is also true for MEBDFI, the FORTRAN implementation of which actually solves index 2 DAEs directly. In addition, the initial condition for the algebraic variable has to be consistent. The problem with using fsolve is that one cannot specify tolerance. Often times one has to solve DAEs at lower tolerances (1e-3 or 1e-4), not 1e-15. In addition, one cannot use compile =true for high efficiency.

The approach in Maple fails for many DAEs of practical interest. In addition, this can give wrong answers. For example, dy/dt = z, y^2+z^2-1=0, with y(0)=0 and z(0)=1 has a meaningful solution only from 0 to Pi/2, Maple’s dsolve/numeric will convert this to dy/dt = z and dz/dt = -y and integrate in time for ever. This is wrong as at t = Pi/2, the system becomes index 2 DAE and there is more than 1 acceptable solution.

We just recently got our paper accepted that helps Maple's dsolve and many ODE solvers in other languages handle DAEs directly. The approach is rather simple, the index 1 DAE is perturbed as dy/dt = f.S ; -mu.diff(g,t) = g. The switch function is a tanh function which is zero till t = tinit (initialization time). Mu is chosen to be a small number. In addition, lhs of the perturbed equation is simplified using approximate initial guesses as Maple cannot handle non-constant Mass matrix. The paper linked below gives below more details.

http://depts.washington.edu/maple/pubs/U_Apprach_FULL_DRAFT.pdf  

Next, I discuss different examples (code attached).

Example 1: Simple DAE (one ODE and one AE), the proposed approach helps solve this system efficiently without knowing the exact initial condition.

Hint: The code is given with a semicolon at the end of the most of the statements for educational purposes for new users. Using semicolon after dsolve numeric in classic worksheet crashes the code (Maplesoft folks couldn’t fix this despite my request).

Example 2:

This is a nickel battery model (one ODE and one AE). This fails with many solvers unless exact IC is given. The proposed approach works well. In particular, stiff=true, implicit=true option is found to be very efficient. The code given in example 1 can be used to solve example 2 or any other DAEs by just entering ODEs, AEs, ICs in proper places.

Example 3:

This is a nonlinear implicit ODE (posted in Mapleprimes earlier by joha, (http://www.mapleprimes.com/questions/203096-Solving-Nonlinear-ODE#answer211682 ). This can be converted to index 1 DAE and solved using the proposed approach.

Example 4:

This example was posted earlier by me in Mapleprimes (http://www.mapleprimes.com/posts/149877-ODEs-PDEs-And-Special-Functions) . Don’t try to solve this in Maple using its dsolve numeric solver for N greater than 5 directly. The proposed approach can handle this well. Tuning the perturbation parameters and using compile =true will help in solving this system more efficiently.

Finally example 1 is solved for different perturbation parameters to show how one can arrive at convergence. Perhaps more sophisticated users and Maplesoft folks can enable this as automatically tuned parameters in dsolve/numeric.

Note:

This post should not be viewed as just being critical on Maple. I have wasted a lot of time assuming that Maple does what it claims in its help file. This post is to bring awareness about the difficulty in dealing with DAEs. It is perfectly fine to not have a DAE solver, but when literature or standard methods are cited/claimed, they should be implemented in proper form. I will forever remain a loyal Maple user as it has enabled me to learn different topics efficiently. Note that Maplesim can solve DAEs, but it is a pain to create a Maplesim model/project just for solving a DAE. Also, events is a pain with Maplesim. The reference to Mapleprimes links are missing in the article, but will be added before the final printing stage. The ability of Maple to find analytical Jacobian helps in developing many robust ODE and DAE solvers and I hope to post my own approaches that will solve more complicated systems.

I strongly encourage testing of the proposed approach and implementation for various educational/research purposes by various users. The proposed approach enables solving lithium-ion and other battery/electrochemical storage models accurately in a robust manner. A disclosure has been filed with the University of Washington to apply for a provisional patent for battery models and Battery Management System for transportation, storage and other applications because of the current commercial interest in this topic (for batteries). In particular, use of this single step avoids intialization issues/(no need to initialize separately) for parameter estimation, state estimation or optimal control of battery models.

Download SolvingDAEsinMaple.mws

I'd like to pay attention to Maple packages for symmetric functions, posets, root systems, and finite Coxeter groups by John Stembridge. That was estimated as well qualified in Mathoverflow.

@Markiyan Hirnyk   I try not to use this package, as I think the results are not reliable enough. Here is the example, where instead of the three real roots it finds only one, despite the hint to look for the three roots:

restart;

DirectSearch:-SolveEquations(100^x=x^100, AllSolutions, solutions=3);

There are many other examples, particularly in discrete optimization in which it returns false results. Here is one example (well-known to you).