- Replies - MaplePrimes

Resolved - at least partially...

Answered: Dr. Venkat Subramanian 381

August 16 2023

0 0

Please see the code attached.
Erik Postma @epostma had helped me in the past to call the MKL solver, so I am just guessing the libraries and trying it. It seems to work. But it may be buggy.

IntelMKL PARDISO can be superior to UMFPACK for very large sparse systems as shown earlier by our group for the paper and codes at https://sites.utexas.edu/maple/phase-field-models-for-li-metal-batteries/

>

restart:

>	Digits:=15;

(1)

>	A:=Matrix(4,4,[[-1,2,0,0],[2,-1,2,0],[0,2,-1,2],[0,0,1,-1]],datatype=float[8],storage=sparse);

Matrix(4, 4, {(1, 1) = -1.0, (1, 2) = 2.0, (1, 3) = 0., (1, 4) = 0., (2, 1) = 2.0, (2, 2) = -1.0, (2, 3) = 2.0, (2, 4) = 0., (3, 1) = 0., (3, 2) = 2.0, (3, 3) = -1.0, (3, 4) = 2.0, (4, 1) = 0., (4, 2) = 0., (4, 3) = 1.0, (4, 4) = -1.0})

(2)

>	b:=Vector(4,[1,0.,0,0],datatype=float[8]):

>

>	sol:=LinearAlgebra:-LinearSolve(A,b,method=SparseDirectMKL);

Error, invalid input: LinearAlgebra:-LinearSolve expects value for keyword parameter method to be of type identical(none,SparseLU,SparseDirect,SparseDirectMKL,SparseIterative,LU,QR,solve,hybrid,Cholesky,subs,modular), but received SparseDirectMKL

>	s1:=LinearAlgebra:-LUDecomposition(A,method=SparseDirectMKL,output=handle);

before P call
after P call

(3)

>

s1();

(4)

>	pSolve := define_external("hw_PardisoSolve", ':-MAPLE', ':-LIB' = "linalg"):

>	x:=Vector(4,datatype=float[8]):

>	pSolve(b,x,s1);

>

x;

Vector(4, {(1) = -5.0, (2) = -2.0, (3) = 4.0, (4) = 4.0})

(5)

>

>	pFactor := define_external("hw_PardisoFactor", ':-MAPLE', ':-LIB' = "linalg"):

>	pCSF := define_external("hw_SpCompressedSparseForm", ':-MAPLE', ':-LIB' = "linalg"):

>

One can also use pFactor, pCSF, pSolve to get C, R, X and then proceed.

>	CB,R,X:=pCSF(A,"r",1,1,false,true,false);

CB, R, X := Vector(5, {(1) = 1, (2) = 3, (3) = 6, (4) = 9, (5) = 11}), Vector(10, {(1) = 1, (2) = 2, (3) = 1, (4) = 2, (5) = 3, (6) = 2, (7) = 3, (8) = 4, (9) = 3, (10) = 4}), Vector(10, {(1) = -1., (2) = 2., (3) = 2., (4) = -1., (5) = 2., (6) = 2., (7) = -1., (8) = 2., (9) = 1., (10) = -1.})

(6)

>	mtype := 11; handle1 := pFactor(CB, R, X, mtype);

before P call
after P call

(7)

>	x:=Vector(4,datatype=float[8]): pSolve(b,x,handle1): x;

Vector(4, {(1) = -5.0, (2) = -2.0, (3) = 4.0, (4) = 4.0})

(8)

Download buglinearsolve.mw

well done...

Answered: Dr. Venkat Subramanian 381

August 02 2023

0 0

Such a great thread, thanks all!

Thanks...

Answered: Dr. Venkat Subramanian 381

June 16 2023

0 0

@Carl Love
Not sure I know the difference. I am just looking for something that will call the external integrator without using evalf command.

See a MATLAB call from Maple and the small amount of RAM used (after the MATLAB opens up the first time). (This calls the vdp1 model). The expected wrapper or procedure should call the external integrator like this for a given n(number of odes), f(procedure to find the slopes), initial condition vector (y0), initial time (to) and final time (tf) with default settings for all other parameters

@AllanWittkopf is trying to see if he can help. I will update this page. My rough test for a sample code was able to run a model with 50KB RAM consumption for Maple's dsolve/numeric external integrator, this involved guessing the entries for the table creation/modification.

>

restart:

>	Digits:=15;

>	currentdir();

$"C:\Program Files\Maple 2023"$

>	with(Matlab);

[AddTranslator, CloseEngine, ExecuteCommand, FromMFile, FromMatlab, GetVariable, IsAssigned, OpenEngine, SetVariable, chol, closelink, defined, det, dimensions, eig, evalM, fft, getvar, inv, lu, ode15s, ode45, openlink, qr, setvar, size, square, transpose]

>	ss:=proc(t0,tf,ic0,T,Y) T,Y := ode45("vdp1", 0 ..tf, ic0); end proc;

>	ic0:=[1,0];

>	CodeTools:-Usage((ss(0,1,ic0,T,Y)));

memory used=0.79MiB, alloc change=0 bytes, cpu time=0ns, real time=4.29s, gc time=0ns

>	T[57];Y[57];

>

Download MATLABcall.mw

Partial answer (not sure) and a potentia...

Answered: Dr. Venkat Subramanian 381

June 14 2023

0 0

Please see the code attached. If anyone can explain the arguments passed then one can possibly write an evalhf-able procedure that will take n (number of odes), f(procedure for the slope), initial time, final time and do all of this with the external integrator and provide the results in evalhf format.
If anyone can show and store the misc and xf called by `dsolve/numeric/SC/IVPrun` this can be done. The goal is to reduce the number of non-evalhf able loops/calls.

(this code is written in Maple 2023).

>

restart:

>	Digits:=15:

>	eqs:=diff(ca(t),t)=-(u+u^2/2)ca(t),diff(cb(t),t)=uca(t)-0.1*cb(t);

>

u:=0.1;

>	infolevel[all]:=10;

>	#soln:=dsolve({eqs,ca(0)=1.0,cb(0)=0.0},type=numeric):

>

#soln(1);

>	#showstat(`dsolve/numeric/SC/IVPrun`);

>	showstat(`dsolve/numeric/SC/hw_int`);

`dsolve/numeric/SC/hw_int` := proc()
local lib;
global `dsolve/numeric/SC/hw_int`;
   1   try
   2       lib := ExternalCalling:-ExternalLibraryName("ode2",'HWFloat');
   3       `dsolve/numeric/SC/hw_int` := define_external('hw_integrate',':-
             MAPLE',(':-LIB') = lib,':-THREAD_SAFE')
       catch :
   4       error
             "'hw_integrate' in external library %1 could not be found/used",
             lib
       end try;
   5   `dsolve/numeric/SC/hw_int`(_passed[1 .. _npassed])
end proc

>

`dsolve/numeric/SC/hw_int` := proc()
local lib;
global `dsolve/numeric/SC/hw_int`;
      try
          lib := ExternalCalling:-ExternalLibraryName("ode2",'HWFloat');
          `dsolve/numeric/SC/hw_int` := define_external('hw_integrate',':-
             MAPLE',(':-LIB') = lib,':-THREAD_SAFE')
       catch :
          error
             "'hw_integrate' in external library %1 could not be found/used",
             lib
       end try;
print(_passed);
      `dsolve/numeric/SC/hw_int`(_passed[1 .. _npassed]);
print("here");
print(_passed);
print(_passed[11]);
end proc;

>	soln:=dsolve({eqs,ca(0)=1.0,cb(0)=0.0},type=numeric,output=Array([seq(i/1.*1,i=0..1)])):

dsolve/numeric: entering dsolve/numeric

DEtools/convertsys: Naming is {ca(t) = Y[1] cb(t) = Y[2] diff(ca(t) t) = YP[1] diff(cb(t) t) = YP[2]}

DEtools/convertsys: converted to first-order system Y'(x) = f(x,Y(x)) namely (with Y' represented by YP)

DEtools/convertsys: correspondence between Y[i] names and original functions:

dsolve/numeric/DAE/rhsproc: Using previously generated rhs proc

dsolve/numeric: the procedure F(x,Y,YP) for computing Y'(x)=f(x,Y(x)) is: proc (N, X, Y, YP) option [Y[1] = ca(t), Y[2] = cb(t)]; YP[1] := -.105000000000000*Y[1]; YP[2] := .1*Y[1]-.1*Y[2]; 0 end proc

dsolve/numeric: initial conditions: x0=0., y0=[1.0, 0.]

dsolve/numeric/SC/preproc: no compilation, _Env_ivpnumeric_ToExternal not set

dsolve/numeric/SC/firststep: Checking ODE for hardware floating point computation

dsolve/numeric/SC/firststep: Initial point OK for hardware floating point computation

dsolve/numeric/SC/firststep: Initial step set to: .0504765875584155

dsolve/numeric/SC/IVPrun: First call, start from initial

dsolve/numeric/SC/IVPrun: Calling external hardware nonstiff integrator with evalhf= true

$Vector[row](2, {(1) = 1., (2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = -.105000000000000, (2) = .100000000000000}), proc (N, X, Y, YP) option `[Y[1] = ca(t), Y[2] = cb(t)]`; YP[1] := (-1)*.105000000000000*Y[1]; YP[2] := .100000000000000*Y[1]+(-1)*.100000000000000*Y[2]; 0 end proc, -1, Vector[row](2, {(1) = .100000000000000, (2) = .100000000000000}), Vector[row](65, {(1) = 2, (2) = 2, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 1, (14) = 2, (15) = 1, (16) = 0, (17) = 0, (18) = 1, (19) = 30000, (20) = 1, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0, (55) = 0, (56) = 0, (57) = 0, (58) = 0, (59) = 10000, (60) = 0, (61) = 1000, (62) = 0, (63) = 0, (64) = -1, (65) = 0}), Vector[row](28, {(1) = 1., (2) = 1.00000000000000*(1/1000000), (3) = 0.504765875584155e-1, (4) = 5.00001000000000*(1/1000000000000000), (5) = 0., (6) = 0.504765875584155e-1, (7) = 0., (8) = 1.00000000000000*(1/1000000), (9) = 0., (10) = 0., (11) = 0., (12) = 0., (13) = 1., (14) = 0., (15) = .499999999999990, (16) = 0., (17) = 1., (18) = 1., (19) = 0., (20) = 0., (21) = 1., (22) = 1., (23) = 0., (24) = 0., (25) = 1.00000000000000*(1/1000000000000000), (26) = 0., (27) = 0., (28) = 0.}), Array(1..6, 0..2, {(1, 1) = .8758457905848228, (1, 2) = .9121381606660928, (2, 0) = .9121381606660928, (2, 1) = 0.8006442103497206e-1, (2, 2) = 1.0962744218225748, (3, 0) = 1.0962744218225748, (3, 1) = .8912691536438646, (3, 2) = 0.9797584316374096e-1, (4, 0) = 0.9797584316374096e-1, (4, 1) = 1.316703053060327, (4, 2) = .8708776141403536, (5, 0) = .8708776141403536, (5, 1) = .11504702599633022, (5, 2) = 1.5371316842980791, (6, 0) = 1.5371316842980791, (6, 1) = .8509526190302513, (6, 2) = .13130657350916902}, datatype = float[8], order = C_order), 0, Array(1..2, 0..2, {(1, 1) = .0, (1, 2) = 1.0, (2, 0) = 1.0, (2, 1) = .0, (2, 2) = 1.0}, datatype = float[8], order = C_order), 0, 0, 0, 0, 0, 0, 0, 0, Matrix(4, 7, {(1, 1) = 0., (1, 2) = .203125000000000, (1, 3) = .304687500000000, (1, 4) = .750000000000000, (1, 5) = .812500000000000, (1, 6) = .406250000000000, (1, 7) = .812500000000000, (2, 1) = 0.637817382812500e-1, (2, 2) = 0., (2, 3) = .279296875000000, (2, 4) = .272378921508789, (2, 5) = -0.968685150146484e-1, (2, 6) = 0.195693969726562e-1, (2, 7) = .538158416748047, (3, 1) = 0.318908691406250e-1, (3, 2) = 0., (3, 3) = -.343750000000000, (3, 4) = -.335235595703125, (3, 5) = .229614257812500, (3, 6) = .417480468750000, (3, 7) = 11.4807128906250, (4, 1) = 0.971052050590515e-1, (4, 2) = 0., (4, 3) = .403503417968750, (4, 4) = 0.202974677085876e-1, (4, 5) = -0.605428218841553e-2, (4, 6) = -0.477004051208496e-1, (4, 7) = .778585672378540}), Matrix(6, 6, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (1, 6) = 1., (2, 1) = .250000000000000, (2, 2) = 0., (2, 3) = 0., (2, 4) = 0., (2, 5) = 0., (2, 6) = 1., (3, 1) = .187500000000000, (3, 2) = .562500000000000, (3, 3) = 0., (3, 4) = 0., (3, 5) = 0., (3, 6) = 2., (4, 1) = .235839843750000, (4, 2) = -.878906250000000, (4, 3) = .890625000000000, (4, 4) = 0., (4, 5) = 0., (4, 6) = .268188476562500, (5, 1) = .127273559570312, (5, 2) = -.500976562500000, (5, 3) = .449218750000000, (5, 4) = -0.128936767578125e-1, (5, 5) = 0., (5, 6) = 0.626220703125000e-1, (6, 1) = -0.927734375000000e-1, (6, 2) = .626220703125000, (6, 3) = -.432617187500000, (6, 4) = .141830444335938, (6, 5) = -0.861053466796875e-1, (6, 6) = .313110351562500}), Vector[row](6, {(1) = 0., (2) = .386000000000000, (3) = .210000000000000, (4) = .630000000000000, (5) = 1., (6) = 1.}), Vector[row](6, {(1) = .250000000000000, (2) = -.104300000000000, (3) = .103500000000000, (4) = -0.362000000000000e-1, (5) = 0., (6) = 0.}), Matrix(6, 5, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (2, 1) = 1.54400000000000, (2, 2) = 0., (2, 3) = 0., (2, 4) = 0., (2, 5) = 0., (3, 1) = .946678528081553, (3, 2) = .255701169898258, (3, 3) = 0., (3, 4) = 0., (3, 5) = 0., (4, 1) = 3.31482518706849, (4, 2) = 2.89612401597212, (4, 3) = .998641913997781, (4, 4) = 0., (4, 5) = 0., (5, 1) = 1.22122450922627, (5, 2) = 6.01913448128775, (5, 3) = 12.5370833293209, (5, 4) = -.687886036105895, (5, 5) = 0., (6, 1) = 1.22122450922627, (6, 2) = 6.01913448128775, (6, 3) = 12.5370833293209, (6, 4) = -.687886036105895, (6, 5) = 1.}), Matrix(6, 5, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (2, 1) = -5.66880000000000, (2, 2) = 0., (2, 3) = 0., (2, 4) = 0., (2, 5) = 0., (3, 1) = -2.43009335683376, (3, 2) = -.206359915708912, (3, 3) = 0., (3, 4) = 0., (3, 5) = 0., (4, 1) = -.107352905814526, (4, 2) = -9.59456225102190, (4, 3) = -20.4702861480962, (4, 4) = 0., (4, 5) = 0., (5, 1) = 7.49644331396861, (5, 2) = -10.2468043146412, (5, 3) = -33.9999035281991, (5, 4) = 11.7089089320616, (5, 5) = 0., (6, 1) = 8.08324679592241, (6, 2) = -7.98113298806279, (6, 3) = -31.5215943287437, (6, 4) = 16.3193054312314, (6, 5) = -6.05881823883405}), Matrix(3, 5, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (2, 1) = 10.1262350834469, (2, 2) = -7.48799587760763, (2, 3) = -34.8009186155574, (2, 4) = -7.99277170756873, (2, 5) = 1.02513772329562, (3, 1) = -.676280339280690, (3, 2) = 6.08771465167861, (3, 3) = 16.4308432089246, (3, 4) = 24.7672251141837, (3, 5) = -6.59438912571678}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Matrix(2, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0.}), Matrix(2, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Matrix(2, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0.}), Matrix(2, 6, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (1, 6) = 0., (2, 1) = 0., (2, 2) = 0., (2, 3) = 0., (2, 4) = 0., (2, 5) = 0., (2, 6) = 0.}), Vector[row](2, {(1) = 0, (2) = 0}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](4, {(1) = 0., (2) = 0., (3) = 0., (4) = 0.}), 0, 0, 0, 0, 0, 0, Vector[row](0, {}), Vector[row](2, {(1) = 0, (2) = 0}), 0$

$Vector[row](2, {(1) = .912138160666093, (2) = 0.800644210349721e-1}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = -0.957745068699397e-1, (2) = 0.832073739631121e-1}), proc (N, X, Y, YP) option `[Y[1] = ca(t), Y[2] = cb(t)]`; YP[1] := (-1)*.105000000000000*Y[1]; YP[2] := .100000000000000*Y[1]+(-1)*.100000000000000*Y[2]; 0 end proc, -1, Vector[row](2, {(1) = .100000000000000, (2) = .100000000000000}), Vector[row](65, {(1) = 2, (2) = 2, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 1, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 1, (14) = 2, (15) = 3, (16) = 0, (17) = 0, (18) = 22, (19) = 30000, (20) = 5, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0, (55) = 0, (56) = 0, (57) = 0, (58) = 0, (59) = 10000, (60) = 0, (61) = 1000, (62) = 0, (63) = 0, (64) = -1, (65) = 0}), Vector[row](28, {(1) = 1., (2) = 1.00000000000000*(1/1000000), (3) = .881714524951008, (4) = 5.00001000000000*(1/1000000000000000), (5) = 0., (6) = 0.504765875584155e-1, (7) = 0., (8) = 1.00000000000000*(1/1000000), (9) = 0., (10) = 0., (11) = 0., (12) = 0., (13) = 1., (14) = 0., (15) = .499999999999990, (16) = 0., (17) = 1., (18) = 1., (19) = 0., (20) = 0., (21) = 1., (22) = 1., (23) = 0., (24) = 0., (25) = 1.00000000000000*(1/1000000000000000), (26) = 0., (27) = 0., (28) = 0.}), Array(1..6, 0..2, {(1, 1) = .8758457905848228, (1, 2) = .9121381606660928, (2, 0) = .9121381606660928, (2, 1) = 0.8006442103497206e-1, (2, 2) = 1.0962744218225748, (3, 0) = 1.0962744218225748, (3, 1) = .8912691536438646, (3, 2) = 0.9797584316374096e-1, (4, 0) = 0.9797584316374096e-1, (4, 1) = 1.316703053060327, (4, 2) = .8708776141403536, (5, 0) = .8708776141403536, (5, 1) = .11504702599633022, (5, 2) = 1.5371316842980791, (6, 0) = 1.5371316842980791, (6, 1) = .8509526190302513, (6, 2) = .13130657350916902}, datatype = float[8], order = C_order), 0, Array(1..2, 0..2, {(1, 1) = .0, (1, 2) = 1.0, (2, 0) = 1.0, (2, 1) = .0, (2, 2) = 1.0}, datatype = float[8], order = C_order), 0, 0, 0, 0, 0, 0, 0, 0, Matrix(4, 7, {(1, 1) = 0., (1, 2) = .203125000000000, (1, 3) = .304687500000000, (1, 4) = .750000000000000, (1, 5) = .812500000000000, (1, 6) = .406250000000000, (1, 7) = .812500000000000, (2, 1) = 0.637817382812500e-1, (2, 2) = 0., (2, 3) = .279296875000000, (2, 4) = .272378921508789, (2, 5) = -0.968685150146484e-1, (2, 6) = 0.195693969726562e-1, (2, 7) = .538158416748047, (3, 1) = 0.318908691406250e-1, (3, 2) = 0., (3, 3) = -.343750000000000, (3, 4) = -.335235595703125, (3, 5) = .229614257812500, (3, 6) = .417480468750000, (3, 7) = 11.4807128906250, (4, 1) = 0.971052050590515e-1, (4, 2) = 0., (4, 3) = .403503417968750, (4, 4) = 0.202974677085876e-1, (4, 5) = -0.605428218841553e-2, (4, 6) = -0.477004051208496e-1, (4, 7) = .778585672378540}), Matrix(6, 6, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (1, 6) = 1., (2, 1) = .250000000000000, (2, 2) = 0., (2, 3) = 0., (2, 4) = 0., (2, 5) = 0., (2, 6) = 1., (3, 1) = .187500000000000, (3, 2) = .562500000000000, (3, 3) = 0., (3, 4) = 0., (3, 5) = 0., (3, 6) = 2., (4, 1) = .235839843750000, (4, 2) = -.878906250000000, (4, 3) = .890625000000000, (4, 4) = 0., (4, 5) = 0., (4, 6) = .268188476562500, (5, 1) = .127273559570312, (5, 2) = -.500976562500000, (5, 3) = .449218750000000, (5, 4) = -0.128936767578125e-1, (5, 5) = 0., (5, 6) = 0.626220703125000e-1, (6, 1) = -0.927734375000000e-1, (6, 2) = .626220703125000, (6, 3) = -.432617187500000, (6, 4) = .141830444335938, (6, 5) = -0.861053466796875e-1, (6, 6) = .313110351562500}), Vector[row](6, {(1) = 0., (2) = .386000000000000, (3) = .210000000000000, (4) = .630000000000000, (5) = 1., (6) = 1.}), Vector[row](6, {(1) = .250000000000000, (2) = -.104300000000000, (3) = .103500000000000, (4) = -0.362000000000000e-1, (5) = 0., (6) = 0.}), Matrix(6, 5, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (2, 1) = 1.54400000000000, (2, 2) = 0., (2, 3) = 0., (2, 4) = 0., (2, 5) = 0., (3, 1) = .946678528081553, (3, 2) = .255701169898258, (3, 3) = 0., (3, 4) = 0., (3, 5) = 0., (4, 1) = 3.31482518706849, (4, 2) = 2.89612401597212, (4, 3) = .998641913997781, (4, 4) = 0., (4, 5) = 0., (5, 1) = 1.22122450922627, (5, 2) = 6.01913448128775, (5, 3) = 12.5370833293209, (5, 4) = -.687886036105895, (5, 5) = 0., (6, 1) = 1.22122450922627, (6, 2) = 6.01913448128775, (6, 3) = 12.5370833293209, (6, 4) = -.687886036105895, (6, 5) = 1.}), Matrix(6, 5, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (2, 1) = -5.66880000000000, (2, 2) = 0., (2, 3) = 0., (2, 4) = 0., (2, 5) = 0., (3, 1) = -2.43009335683376, (3, 2) = -.206359915708912, (3, 3) = 0., (3, 4) = 0., (3, 5) = 0., (4, 1) = -.107352905814526, (4, 2) = -9.59456225102190, (4, 3) = -20.4702861480962, (4, 4) = 0., (4, 5) = 0., (5, 1) = 7.49644331396861, (5, 2) = -10.2468043146412, (5, 3) = -33.9999035281991, (5, 4) = 11.7089089320616, (5, 5) = 0., (6, 1) = 8.08324679592241, (6, 2) = -7.98113298806279, (6, 3) = -31.5215943287437, (6, 4) = 16.3193054312314, (6, 5) = -6.05881823883405}), Matrix(3, 5, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (2, 1) = 10.1262350834469, (2, 2) = -7.48799587760763, (2, 3) = -34.8009186155574, (2, 4) = -7.99277170756873, (2, 5) = 1.02513772329562, (3, 1) = -.676280339280690, (3, 2) = 6.08771465167861, (3, 3) = 16.4308432089246, (3, 4) = 24.7672251141837, (3, 5) = -6.59438912571678}), Vector[row](2, {(1) = .862852234180585, (2) = .121646667791134}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = -0.873057666308174e-1, (2) = 0.684701246507017e-1}), Vector[row](2, {(1) = 0., (2) = 0.}), Matrix(2, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0.}), Matrix(2, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Matrix(2, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0.}), Matrix(2, 6, {(1, 1) = -0.957745068699397e-1, (1, 2) = -0.935578053069745e-1, (1, 3) = -0.925071732935370e-1, (1, 4) = -0.879217450216900e-1, (1, 5) = -0.872861411194737e-1, (1, 6) = -0.914468108408395e-1, (2, 1) = 0.832073739631121e-1, (2, 2) = 0.792621008622733e-1, (2, 3) = 0.774422703214511e-1, (2, 4) = 0.695044462636171e-1, (2, 5) = 0.684014621613813e-1, (2, 6) = 0.755991753809135e-1}), Vector[row](2, {(1) = 0, (2) = 0}), Vector[row](2, {(1) = .912138160666093, (2) = 0.800644210349721e-1}), Vector[row](2, {(1) = .831483491722071, (2) = .146782245215054}), Vector[row](2, {(1) = 1.66780067534233*(1/1000000000), (2) = 6.93290889974829*(1/1000000000)}), Vector[row](2, {(1) = .912138160666093, (2) = 0.800644210349721e-1}), Vector[row](2, {(1) = -0.903900345181292e-1, (2) = 0.737580972195178e-1}), Vector[row](4, {(1) = 0., (2) = 0., (3) = 0., (4) = 0.}), 0, 0, 0, 0, 0, 0, Vector[row](0, {}), Vector[row](2, {(1) = 0, (2) = 0}), 0$

Matrix(2, 3, {(1, 1) = 0., (1, 2) = 1.000000000000000, (1, 3) = 0., (2, 1) = 1.000000000000000, (2, 2) = .9003245213933421, (2, 3) = 0.9025791441096123e-1})

dsolve/numeric/SC/IVPrun: Completed external hardware nonstiff integrator

dsolve/numeric/SC/IVPrun: Time for integration= 0.

>

What we need to do is, write a procedure that will take inputs as the number of odes, n, procedure f, initial conditions for Y vector, initial time, final time, and update the arguments and provide the final Y value from the arguments in evalhf form by direcly calling `dsolve/numeric/SC/hw_int`(_passed[1 .. _npassed]);

Download learningdsolvenumeric.mw

bump up...

Answered: Dr. Venkat Subramanian 381

June 14 2023

0 0

Trying to bump up this thread to get the attention of @acer, @tomleslie, and others.
Hopefully, this example convinces everyone that just choosing a different optimization package (global optimization, stochastic, genetic algorithm, etc) may not help in achieving the optimal results if the numerical discretization (in this case nvar for the number time horizons) is not providing sufficient precision. We would need nvar>=256 for this problem and approach.

If anyone has access to old Maple versions, calling the external dsolve/numeric (Without the print statements) and getting an output array as shown in the link below and example towards the end for the NAG subroutine should help.

The link is at
https://www.maplesoft.com/support/helpJP/Maple/view.aspx?path=NAG/d02cjc

Though Maple does not support NAG toolbox anymore, if anyone has access to Maple 10 or earlier and if you are able to run the example on the NAG help page, please let me know.

Thanks

GlobalSolve...

Answered: Dr. Venkat Subramanian 381

June 05 2023

0 0

@sursumCorda

This is based on my limited testing of combining Maple's code with its global optimization toolbox.

Currently, you pay a separate license for the same.

Currently, this is offered by Nimbus which provides Genetic Algorithms.
A decade or so earlier, Maple provided the Global Optimization toolbox from Pintur. That had method=reducedgradient which was very robust and fast.
Hope this helps.

Thanks and a general comment as a profes...

Answered: Dr. Venkat Subramanian 381

June 05 2023

0 0

@tomleslie

I like DirectSearch and I have nothing against it or genetic algorithms (GA) which are heavily used in the literature.

This discussion and analysis (from @acer and the solution from him) is the best and most efficient way to solve this model.

The next biggest step can be wrapping dsolve/numeric (with or without compile=true) and making the NLPSolve work with the same in evalhf form for computational efficiency.

The last 5-10 years have seen an exponential increase in open access/open software codes in Python, Julia, etc. Many times I have wondered why students/researchers used GA for well-defined models (with convexity). Now, I partially understand the reason. Many times wrappers are written for optimizers or different algorithms and then called from user-friendly software. But the wrapper is probably poorly written so that the gradient information is not properly passed with a sufficient number of options (tolerances in this case).

PS, note that only gradient-based methods have theoretical backing for locally optimal results for optimization and I am weary of using methods that don't have a good theoretical basis, though they might work very well (but they will be always slower for problems that can be solved with gradient based methods).

Thanks...

Answered: Dr. Venkat Subramanian 381

June 02 2023

0 0

@acer

I think any opimizer or nonlinearsolver should take the required tolerance as input and return the current-best objective, iterate.
@acer, I can confirm that this approach works with nonlinear constraints + dsolve/numeric/parameters + sqp. thanks again

Thanks, this is the best answer...

Answered: Dr. Venkat Subramanian 381

June 01 2023

0 0

@acer

There is no point in doing optimization to 1e-16 accuracy when we are trying to increase nvar(number of stages), solving ODEs/DAEs to 1e-6 accuracy.

I like your answer and it works.

Note that for some weird reason, for this problem stiff=true in the call to dsolve/numeric works without tweaking anything in the dsolve command or optimization command.

Knowing where it fails can help. For example, ideally NLPSolve should have given the answer that "objective was solved to 1e-6 precision, but not beyond that", instead of just returning "initial point satisfies first-order optimality condition".

fsolve has the same issue as well, and if I am not wrong, we can't force it to solve to a prescribed precision.

@acer, as again, thanks lot. I believe that you can fix almost any issues!

Found a solution for now- still iffy...

Answered: Dr. Venkat Subramanian 381

June 01 2023

0 0

Adding abstol=1e-12, reltol=1e-12 to dsolve/numeric calls solves the problem for this model. (Thanks to Allan Wittkopf)

He had also suggested using Digits:=18 inside ss (I am not a big fan of this)

With default tolerances use of, maxstep = 0.05 (assuming that your maximum time horizon is 0.1 with nvar = 10)seems to work as well.

fdiff or the procedures used by NLPSolve should be made more robust.

Thanks...

Answered: Dr. Venkat Subramanian 381

May 31 2023

0 0

@dharr

Your code is not using dsolve numeric in parametric form. I will be posting the use of parametric form working with NLPSolve (thanks to Allan Wittkopf).

As I said, we don't want analytical solutions.

Also, one parameter optimization is different from multiple parameters (scalar vs vector).

To be clear, I am not looking for an answer for the model. My goal was to use NLPSolve with dsolve/numeric/parameters with sqp.

this is not parametric dsolve/numeric...

Answered: Dr. Venkat Subramanian 381

May 31 2023

0 0

@dharr

thanks

Note how the code works when dsolve nume...

Answered: Dr. Venkat Subramanian 381

May 31 2023

0 0

When dsolve numeric is replaced with the RK2 method, the code works. Typically, dsolve numeric will be faster than writing for loops to do numerical integration of ODEs.

This is why I think that the issue is with dsolve/numeric.

>

restart:

>

Test code written by Dr. Venkat Subramanian at UT Austin, 05/31/2023. This code uses CVP approach (piecwise constant) to perform optimal control. NLPSolve combination with RK2 approach to integrate ODEs (constant step size) works well for this problem.

>	Digits:=15;

>

>	eqodes:=[diff(ca(t),t)=-(u+u^2/2)1.0ca(t),diff(cb(t),t)=1.0uca(t)-0.1*cb(t)];

[diff(ca(t), t) = -1.0*(u+(1/2)*u^2)*ca(t), diff(cb(t), t) = 1.0*u*ca(t)-.1*cb(t)]

>	soln:=dsolve({op(eqodes),ca(0)=alpha,cb(0)=beta},type=numeric,'parameters'=[alpha,beta,u],compile=true,savebinary=true):

>

RK2:=proc(alpha,beta,NN,u,dt)#this procedure can be made efficient in vector form
local j,c1mid,c2mid,c10,c20,c1,c2;
c10:=alpha;c20:=beta;
for j from 1 to NN do
  c1mid:=c10+dt/NN*(-(u+u^2/2)*c10):
  c2mid:=c20+dt/NN*(u*c10)-0.1*dt/NN*c20:
  c1:=c10/2+c1mid/2+dt/NN/2*(-(u+u^2/2)*c1mid):
  c2:=c20/2+c2mid/2+dt/NN/2*(u*c1mid)-0.1*dt/NN/2*c2mid:
  c10:=c1:c20:=c2:
  od:
[c10,c20];
end proc:

>	soln('parameters'=[1,0,0.1]);soln(0.1);

>	RK2(1,0,20,0.1,0.1);#reasonable accuracy with just 20 time steps of RK2

>

ss0:=proc(x)
interface(warnlevel=0):
#if type(x[1],numeric)
if type(x,Vector)
then

local z1,n1,i,c10,c20,dt,u;
global soln,nvar;
dt:=evalf(1.0/nvar):
c10:=1.0:c20:=0.0:
for i from 1 to nvar do
u:=x[i]:
soln('parameters'=[c10,c20,u]):
z1:=soln(dt):
c10:=subs(z1,ca(t)):c20:=subs(z1,cb(t)):
od:
-c20;
else 'procname'(args):

end if:

end proc:

>

ss:=proc(x)#based on RK2
interface(warnlevel=0):
#if type(x[1],numeric)
if type(x,Vector)
then

local z1,n1,i,c10,c20,dt,u,NN;
global soln,nvar;
dt:=evalf(1.0/nvar):
NN:=20;
c10:=1.0:c20:=0.0:
for i from 1 to nvar do
u:=x[i]:
z1:=RK2(c10,c20,NN,u,dt):
c10:=z1[1]:c20:=z1[2]:
od:
-c20;
else 'procname'(args):

end if:

end proc:

>	nvar:=3;#code works for nvar:=3, but not for nvar:=2

>	ic0:=Vector(nvar,[seq(0.1,i=1..nvar)],datatype=float[8]);

Vector(3, {(1) = .1000000000000000, (2) = .1000000000000000, (3) = .1000000000000000})

>

ss(ic0);

>	bl := Vector(nvar,[seq(0.,i=1..nvar)]);bu := Vector(nvar,[seq(5.,i=1..nvar)]);

Vector(3, {(1) = 0., (2) = 0., (3) = 0.})

Vector[column](%id = 36893490450763924652)

>	infolevel[all]:=1;

>	Optimization:-NLPSolve(nvar,ss,[],initialpoint=ic0,[bl,bu]);

NLPSolve: calling NLP solver

NLPSolve: using method=sqp

NLPSolve: number of problem variables 3

NLPSolve: number of nonlinear inequality constraints 0

NLPSolve: number of nonlinear equality constraints 0

NLPSolve: number of general linear constraints 0

[-.531367180853919, Vector(3, {(1) = .8116892515974860, (2) = 1.189817750741097, (3) = 2.426110609252688})]

>	nvar:=2;#code works for nvar:=3, but not for nvar:=2

>	ic0:=Vector(nvar,[seq(0.1,i=1..nvar)],datatype=float[8]);

Vector(2, {(1) = .1000000000000000, (2) = .1000000000000000})

>

ss(ic0);

>	bl := Vector(nvar,[seq(0.,i=1..nvar)]);bu := Vector(nvar,[seq(5.,i=1..nvar)]);

Vector(2, {(1) = 0., (2) = 0.})

Vector[column](%id = 36893490450763939460)

>	infolevel[all]:=1;

>	Optimization:-NLPSolve(nvar,ss,[],initialpoint=ic0,[bl,bu]);

NLPSolve: calling NLP solver

NLPSolve: using method=sqp

NLPSolve: number of problem variables 2

NLPSolve: number of nonlinear inequality constraints 0

NLPSolve: number of nonlinear equality constraints 0

NLPSolve: number of general linear constraints 0

[-.523759339206159, Vector(2, {(1) = .9103276760025985, (2) = 1.925973211795306})]

>

>	nvar:=5;#code works for nvar:=3, but not for nvar:=2

>	ic0:=Vector(nvar,[seq(0.1,i=1..nvar)],datatype=float[8]);

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

>

ss(ic0);

>	bl := Vector(nvar,[seq(0.,i=1..nvar)]);bu := Vector(nvar,[seq(5.,i=1..nvar)]);

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

Vector[column](%id = 36893490450873412892)

>	infolevel[all]:=1;

>	Optimization:-NLPSolve(nvar,ss,[],initialpoint=ic0,[bl,bu]);

NLPSolve: calling NLP solver

NLPSolve: using method=sqp

NLPSolve: number of problem variables 5

NLPSolve: number of nonlinear inequality constraints 0

NLPSolve: number of nonlinear equality constraints 0

NLPSolve: number of general linear constraints 0

[-.537292701143927, Vector(5, {(1) = .7437459723532210, (2) = .9015469842007691, (3) = 1.149056357665623, (4) = 1.629519387333676, (5) = 3.220718773451731})]

>	infolevel[all]:=0:

>	#Note that results from nvar =2 are good initial guesses for nvar =4, etc, but it is not implemented in this code.

>	# objective is approaching convergence with 16 stages for control variable.

>	for j from 1 to 4 do nvar:=2^j: ic0:=Vector(nvar,[seq(0.1,i=1..nvar)],datatype=float[8]); bl := Vector(nvar,[seq(0.,i=1..nvar)]);bu := Vector(nvar,[seq(5.,i=1..nvar)]); obj[j]:=Optimization:-NLPSolve(nvar,ss,[],initialpoint=ic0,[bl,bu]): print(2^j,obj[j][1]): gc(); od:

>	[seq(obj[4][2][i],i=1..16)];#Control values

>

Download testRk2.mw

request and apologies...

Answered: Dr. Venkat Subramanian 381

May 31 2023

0 0

@dharr

I am a long-time Maple user, and I would appreciate it if you can provide any code/answer without the shortcuts/etc as much as possible. Also, for any procedures please don't add any inputs/variables/descriptions/variable types/ etc unless they are critically needed to make the code work for that purpose.

One of the charms of Maple is its symbolic capability, and keeping the codes/commands minimal and getting the results we seek are the best options for people like me who still miss the classic worksheet and old Maple GUIs.

Thanks

Venkat

thanks and sorry...

Answered: Dr. Venkat Subramanian 381

May 31 2023

0 0

@dharr

As stated early, I am only interested in the sqp and gradient-based method. Nonlinear simplex cannot handle bounds.
My interest lies in sqp as I have to be able to able to change nvar as input to a procedure to predict the objective, and possibly nonlinear constraints as an objective

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

Venkat Subramanian

381 Reputation

13 Badges

MaplePrimes Activity

These are replies submitted by

Resolved - at least partially...

well done...

Thanks...

Partial answer (not sure) and a potentia...

bump up...

GlobalSolve...

Thanks and a general comment as a profes...

Thanks...

Thanks, this is the best answer...

Found a solution for now- still iffy...

Thanks...

this is not parametric dsolve/numeric...

Note how the code works when dsolve nume...

request and apologies...

thanks and sorry...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again