Advanced Search

Venkat Subramanian

386 Reputation

13 Badges

15 years, 242 days
Contact 119368

MaplePrimes Activity


These are replies submitted by

bump up...

Dr. Venkat Subramanian 386
April 27 2022
0 0

@acer

bump up. Any luck with this? Is it possible to store the compiled procedures in a library or a code so that they can be used for future uses?

thanks

VS

of course Linux is ok to start with...

Dr. Venkat Subramanian 386
April 23 2022
0 0

@ acer

I am ok if you start with Linux. We can always create codes and run them later in linux. But we have to provide a way for the users to run the code with pardiso solver (and umfpack solver) in a smooth manner so that the codes can be used without too much frustration.

After 10^4*10^4 iterated solvers may be faster, but direct solvers are very robust, and may be needed for these problems.

Windows please...

Dr. Venkat Subramanian 386
April 23 2022
0 0

@acer 

Thanks. I was able to run your code. The wrapper you wrote was very helpful for me earlier for Rosenbrock methods (LU decomposition is done once with multiple solves). 
If you run pardiso codes for N =160, 200 or more, you will see more gain from pardiso compared to the umfpack code.

I ran the standard code (umfpack code) in Maple 2020 (it works in 2021 or 2022 as well).

The pardiso code was run in Maple 2021.1 in windows version.  

I prefer the Windows version, as I start building my codes using the classic worksheet Maple (in Maple 14 or so) and then move to mw format. 

For the SLU wrapper, why do we need SLU can't we directly use 

Anz,NumericA := HWcall[1](nrowsA,ncolsA,A);
  res := HWcall[2](nrowsA,m,Anz,NumericA,localx);

after defing HWcall[1], call[2], etc?

 

thanks a lot

VS


 

2 of the 3 suggestions...

Dr. Venkat Subramanian 386
April 22 2022
0 0

@ acer,

I am not sure. I believe I may have already implemented the first 2 of your 3 suggestions in the pardiso code.

 

Excellent...

Dr. Venkat Subramanian 386
April 22 2022
0 0

@acer 

Awesome. Please do so. Thanks. Phase-field models might require 500 or more nodes in each dimension for sufficient precision. Please find enclosed a pardiso based linear solver code (in which I had optimized the calls to the sparse matrix storage using just the non-zero entries). This code will not directly run in Maple, unless you get the MKL specific dll setup properly, but you can see the structure of the code here. I strongly believe UMFPACK can be used similarly (based on the instructions given by UMFPACK developers).

In the paper, we also introduce new RK methods called approximate RK3a in which linearsolver is called only once during the first stage, with interpolations at the other stages.

The Pardiso based Maple code is attached below. We were able to run 1000x1000 or more in an hour or so. This code also defines all the necessary procedures and uses the same compiled procedures to simulate the results for a different number of node points.

Most of these procedures can be preloaded as a separate code (or module) and directly called with N and design parameters in the main Maple window. But I want to give the codes with a minimum number of commands and full transparency.
 

> 

restart:
> 

 

> 

 

04/21/2022. Code written by Dr. Venkat Subramanian and MAPLE group at UT. The code is provided as open-access with no restrictions.

This code is very similar to the standard 2D Battphase code. The only difference is the Pardiso linear solver is called after changing the dll files. Both the clock time and total computation time are reported. CPU reported in Table 3 at the end of the code for  160 node points in each direction. The entire Jacobian is not used. Only the non-zero sparse entries are identified (only once) and updated with time.

> 

t11:=time[real]():t12:=time():
> 

gc();
> 

Digits:=15;

15

 (1)
> 

 

> 

interface(prettyprint=2);

2

 (2)
> 

EFAdd:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float[8],N::integer,M::integer,e1::float[8])
local i::integer,j::integer;
for i from 2 to N-1 do for j from 2 to M-1 do Ymid[i,j]:=max(e1,Y00[i,j]-dt*F0[i,j]*Phi0[i+(j-1)*N]):od:od:
for i from 1 to N do Ymid[i,1]:=Ymid[i,2]:Ymid[i,M]:=Ymid[i,M-1]: od:
for j from 1 to M do Ymid[1,j]:=Ymid[2,j]:Ymid[N,j]:=Ymid[N-1,j]: od:
end proc:
> 

EFAdd:=Compiler:-Compile(EFAdd):
> 

EFAdd2:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float[8],N::integer,M::integer,e1::float[8])
local i::integer,j::integer;
for i from 2 to N-1 do for j from 2 to M-1 do
Ymid[i,j]:=max(e1,Y00[i,j]*3/4.+Ymid[i,j]/4.-dt/4.*F0[i,j]*Phi0[i+(j-1)*N]):od:od:
for i from 1 to N do Ymid[i,1]:=Ymid[i,2]:Ymid[i,M]:=Ymid[i,M-1]:od:
for j from 1 to M do Ymid[1,j]:=Ymid[2,j]:Ymid[N,j]:=Ymid[N-1,j]:od:
end proc:
> 

EFAdd2:=Compiler:-Compile(EFAdd2):
> 

EFAdd3:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float[8],N::integer,M::integer,e1::float[8])
local i::integer,j::integer;
for i from 2 to N-1 do for j from 2 to M-1 do
Y00[i,j]:=max(e1,Y00[i,j]*1/3.+Ymid[i,j]*2/3.-dt*2/3.*F0[i,j]*Phi0[i+(j-1)*N]):od:od:
for i from 1 to N do Y00[i,1]:=Y00[i,2]:Y00[i,M]:=Y00[i,M-1]:od:
for j from 1 to M do Y00[1,j]:=Y00[2,j]: Y00[N,j]:=Y00[N-1,j]:od:
end proc:
> 

EFAdd3:=Compiler:-Compile(EFAdd3):
> 

YdatStore:=proc(Y00::Matrix(datatype=float[8]),Ydat::Array(datatype=float[8]),N::integer,M::integer,jj::integer)
local i::integer,j::integer;
for i from 2 to N-1 do for j from 2 to M-1 do
Ydat[jj,i,j]:=Y00[i,j]:od:od:
end proc:
> 

YdatStore:=Compiler:-Compile(YdatStore):
> 

phiaveAdd:=proc(Y00::Matrix(datatype=float[8]),N::integer,M::integer,Ny::integer)
local i::integer,j::integer,phiave::float;
phiave:=0.0:
for i from 2 to N-1 do for j from 2 to M-1 do phiave:=phiave+Y00[i,j]:od:od:
phiave/(N-2)/(M-2);
end proc:
> 

phiaveAdd:=Compiler:-Compile(phiaveAdd):
> 

PhiAdd:=proc(N::integer,Phi0::Vector(datatype=float[8]),db::Vector(datatype=float[8]))
local i::integer;
for i from 1 to N do Phi0[i]:=Phi0[i]+db[i]:od:
end proc:
> 

PhiAdd:=Compiler:-Compile(PhiAdd):
> 

PhiSwitch:=proc(N::integer,Phi0::Vector(datatype=float[8]),Phitemp::Vector(datatype=float[8]))
local i::integer;
for i from 1 to N do Phitemp[i]:=Phi0[i]:od:
end proc:
> 

PhiSwitch:=Compiler:-Compile(PhiSwitch):
> 

MidPred:=proc(N::integer,Phi0::Vector(datatype=float[8]),Phitemp::Vector(datatype=float[8]),Phimid::Vector(datatype=float[8]))
local i::integer;
for i from 1 to N do Phimid[i]:=Phitemp[i]/2.0+Phi0[i]/2.0:od:
end proc:
> 

MidPred:=Compiler:-Compile(MidPred):
> 

y0proc:=proc(NN::integer,Y00::Matrix(datatype=float[8]),beta::float[8],e1::float[8])
local i,j,xx,yy,N,h,w,MM,M,Ny,rf,f,ff,rr;
> 

N:=NN+2:h:=1.0/NN:w:=h*beta:
> 

MM:=NN;M:=MM+2;Ny:=MM;
> 

for i from 2 to N-1 do for j from 2 to M-1 do
xx:=-0+(i-1/2-1)*h:yy:=-0+(j-1/2-1)*h:
rr:=xx^2+yy^2;
#if j<4 and i <N/5+1 then Y00[i,j]:=1e-9 else Y00[i,j]:=1.0;end:
#Y00[i,j]:=1.0:
#if rr <9/100. then Y00[i,j]:=1e-12; else Y00[i,j]:=1.0:end:
#if rr=9/100. then Y00[i,j]:=0.5:end:
Y00[i,j]:=max(e1,0.5+0.5*tanh((sqrt(rr)-0.3)/w/sqrt(2.0))):
#if abs(xx-0.5)>=0.25 and yy<=0.5 then Y00[i,j]:=1e-9 else Y00[i,j]:=1.0:end:
#if yy<=0.1 then Y00[i,j]:=1e-9 :end:
od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:
> 

 

> 

y0proc0:=proc(NN,Y00)
local i,j,xx,yy,N,h,w,MM,M,Ny,rf,f,ff,rr;
> 

N:=NN+2:h:=1.0/NN:w:=h/2.0:
> 

MM:=NN;M:=MM+2;Ny:=MM;
> 

for i from 2 to N-1 do for j from 2 to M-1 do
xx:=-0+(i-1/2-1)*h:yy:=-0+(j-1/2-1)*h:
rr:=xx^2+yy^2;
#if j<4 and i <N/5+1 then Y00[i,j]:=1e-9 else Y00[i,j]:=1.0;end:
#Y00[i,j]:=1.0:
#if rr <9/100. then Y00[i,j]:=1e-12; else Y00[i,j]:=1.0:end:
#if rr=9/100. then Y00[i,j]:=0.5:end:
#Y00[i,j]:=max(1e-12,0.5+0.5*tanh((sqrt(rr)-0.3)/w/sqrt(2.0))):
if abs(xx-0.5)>=0.25 and yy<=0.5 then Y00[i,j]:=1e-9 else Y00[i,j]:=1.0:end:
if yy<=0.1 then Y00[i,j]:=1e-9 :end:
od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:
> 

y0proc:=Compiler:-Compile(y0proc):
> 

printlevel:=2;

2

 (3)
> 

 

> 

Eqs11:=proc(N::integer,M::integer,Y0::Matrix(datatype=float[8]),F::Matrix(datatype=float[8]),y::Vector(datatype=float[8]),v0::float,ff::Vector(datatype=float[8]))
local i::integer,j::integer,i1::integer,h::float[8];
h:=1.0/(N-2):
j:=1: for i from 1 to N do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1+N]:
od:

for j from 2 to M-1 do
i:=1:
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1+1]:
for i from 2 to N-1 do
i1:=i+(j-1)*N:
ff[i1]:=
(Y0[i,j]+Y0[i,j+1])*(y[i1+N]-y[i1])
-(Y0[i,j]+Y0[i,j-1])*(y[i1]-y[i1-N])
+(Y0[i+1,j]+Y0[i,j])*(y[i1+1]-y[i1])
-(Y0[i,j]+Y0[i-1,j])*(y[i1]-y[i1-1])
-y[i1]*h*(1e-24+(Y0[i+1,j]-Y0[i-1,j])^2+(Y0[i,j+1]-Y0[i,j-1])^2)^(1/2):
od:
i:=N:
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1-1]+0*h:
od:

j:=M:
for i from 1 to N do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1-N]+v0*h:
od:

end proc:
> 

Eqs11:=Compiler:-Compile(Eqs11):
> 

 

> 

XX:=proc(N::integer,M::integer,Y0::Matrix(datatype=float[8]),y::Vector(datatype=float[8]),v0::float,x0::Vector(datatype=float[8]))
local i::integer,j::integer,i1::integer,h::float[8],count::integer;
h:=1.0/(N-2):
count:=0:
j:=1:
for i from 1 to N do
i1:=i+(j-1)*N:
#ff[i1]:=-y[i1]+y[i1+N]:
count:=count+1:
x0[count]:=1.0:
count:=count+1:
x0[count]:=-1.0:
#j00[i1,i1]:=1.0:j00[i1,i1+N]:=-1.00:
od:

for j from 2 to M-1 do
i:=1:
i1:=i+(j-1)*N:
count:=count+1:
x0[count]:=1.0:
count:=count+1:
x0[count]:=-1.0:
#j00[i1,i1]:=1.0:j00[i1,i1+1]:=-1.00:
for i from 2 to N-1 do
i1:=i+(j-1)*N:
count:=count+1:
x0[count]:=-Y0[i,j-1]-Y0[i,j]:
count:=count+1:
x0[count]:=-Y0[i-1,j]-Y0[i,j]:
count:=count+1:
x0[count]:=4*Y0[i,j]+Y0[i,j+1]+Y0[i,j-1]+Y0[i+1,j]+Y0[i-1,j]+h*(1e-24+(Y0[i+1,j]-Y0[i-1,j])^2+(Y0[i,j+1]-Y0[i,j-1])^2)^(1/2):
count:=count+1:
x0[count]:=-Y0[i+1,j]-Y0[i,j]:
count:=count+1:
x0[count]:=-Y0[i,j+1]-Y0[i,j]:
#j00[i1,i1]:=4*Y0[i,j]+Y0[i,j+1]+Y0[i,j-1]+Y0[i+1,j]+Y0[i-1,j]+h*(1e-24+(Y0[i+1,j]-Y0[i-1,j])^2+(Y0[i,j+1]-Y0[i,j-1])^2)^(1/2):
#j00[i1,i1+1]:=-Y0[i+1,j]-Y0[i,j]:j00[i1,i1-1]:=-Y0[i-1,j]-Y0[i,j]:
#j00[i1,i1+N]:=-Y0[i,j+1]-Y0[i,j]:j00[i1,i1-N]:=-Y0[i,j-1]-Y0[i,j]:
od:
i:=N:
count:=count+1:
x0[count]:=-1.0:
count:=count+1:
x0[count]:=1.0:
i1:=i+(j-1)*N:
#j00[i1,i1]:=1.0:j00[i1,i1-1]:=-1.00:
od:

j:=M:
for i from 1 to N do
i1:=i+(j-1)*N:
count:=count+1:
x0[count]:=-1.0:
count:=count+1:
x0[count]:=1.0:
#j00[i1,i1]:=1.0:j00[i1,i1-N]:=-1.00:
#ff[i1]:=-y[i1]+y[i1-N]+v0*h:
od:

end proc:
> 

XX:=Compiler:-Compile(XX):
> 

CBproc:=proc(N::integer,M::integer,CB0::Vector(datatype=integer[4]))
local i::integer,j::integer;
CB0[1]:=1:
> 

for i from 1 to N do CB0[i+1]:=CB0[i]+2:od:
> 

for j from 2 to M-1 do
i:=1:CB0[i+(j-1)*N+1]:=CB0[i+(j-1)*N]+2:
for i from 2 to N-1 do
CB0[i+(j-1)*N+1]:=CB0[i+(j-1)*N]+5:od:
i:=N:CB0[i+(j-1)*N+1]:=CB0[i+(j-1)*N]+2:od:

> 

for i from 1 to N do CB0[N*M-N+1+i]:=CB0[N*M-N+i]+2:od:
end proc:
> 

CBproc:=Compiler:-Compile(CBproc):
> 

Rproc:=proc(N::integer,M::integer,R0::Vector(datatype=integer[4]))
local i::integer,j::integer,i1::integer,count::integer;
count:=0:
j:=1:
for i from 1 to N do
i1:=i+(j-1)*N:
#ff[i1]:=-y[i1]+y[i1+N]:
count:=count+1:
R0[count]:=i1:
count:=count+1:
R0[count]:=i1+N:
#j00[i1,i1]:=1.0:j00[i1,i1+N]:=-1.00:
od:

for j from 2 to M-1 do
i:=1:
i1:=i+(j-1)*N:
count:=count+1:
R0[count]:=i1:
count:=count+1:
R0[count]:=i1+1:
#j00[i1,i1]:=1.0:j00[i1,i1+1]:=-1.00:
for i from 2 to N-1 do
i1:=i+(j-1)*N:
count:=count+1:
R0[count]:=i1-N:
count:=count+1:
R0[count]:=i1-1:
count:=count+1:
R0[count]:=i1:
count:=count+1:
R0[count]:=i1+1:
count:=count+1:
R0[count]:=i1+N:
#j00[i1,i1]:=4*Y0[i,j]+Y0[i,j+1]+Y0[i,j-1]+Y0[i+1,j]+Y0[i-1,j]+h*(1e-24+(Y0[i+1,j]-Y0[i-1,j])^2+(Y0[i,j+1]-Y0[i,j-1])^2)^(1/2):
#j00[i1,i1+1]:=-Y0[i+1,j]-Y0[i,j]:j00[i1,i1-1]:=-Y0[i-1,j]-Y0[i,j]:
#j00[i1,i1+N]:=-Y0[i,j+1]-Y0[i,j]:j00[i1,i1-N]:=-Y0[i,j-1]-Y0[i,j]:
od:
i:=N:
i1:=i+(j-1)*N:
count:=count+1:
R0[count]:=i1-1:
count:=count+1:
R0[count]:=i1:
i1:=i+(j-1)*N:
#j00[i1,i1]:=1.0:j00[i1,i1-1]:=-1.00:
od:

j:=M:
for i from 1 to N do
i1:=i+(j-1)*N:
count:=count+1:
R0[count]:=i1-N:
count:=count+1:
R0[count]:=i1:
#j00[i1,i1]:=1.0:j00[i1,i1-N]:=-1.00:
#ff[i1]:=-y[i1]+y[i1-N]+v0*h:
od:

end proc:
> 

Rproc:=Compiler:-Compile(Rproc):
> 

 

> 

 

> 

ENO2:=proc(Y00::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,N::integer,M::integer,Nx::integer,v0::float)
local i::integer,j::integer,h::float[8],nx::float[8],vel::float[8],vx::float[8],vy::float[8],phix::float[8],phiy::float[8],phiave::float[8],jj::integer,phixb::float[8],phixf::float[8],phixb2::float[8],phixf2::float[8],vxb::float[8],phiyb::float[8],phiyf::float[8],vxf::float[8],phiyb2::float[8],phiyf2::float[8],vyb::float[8],vyf::float[8],uf::float[8],ub::float[8],vf::float[8],vb::float[8],tt::float[8],vv0::float[8],sd::float[8],sdf::float[8],sdb::float[8],sdx::float[8],sdy::float[8],sdxb::float[8],s1x::float[8],sdxf::float[8],sdyb::float[8],sdyf::float[8],s1y::float[8],vx1::float[8],vx2::float[8],vy1::float[8],vy2::float[8],alpha::float[8],beta::float[8];


h:=1.0/(N-2):
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
for i from 2 to N-1 do for j from 2 to M-1 do
vx:=0.0:vy:=0.0:phix:=0.0:phiy:=0.0:

sdx:=(Y00[i+1,j]-2*Y00[i,j]+Y00[i-1,j])/h:sdy:=(Y00[i,j+1]-2*Y00[i,j]+Y00[i,j-1])/h:
vxb:=0:vxf:=0:vyb:=0:vyf:=0:
if i = 2 then sdxb:=(Y00[i,j]-2*Y00[i-1,j]+Y00[i-1,j])/h:else sdxb:=(Y00[i,j]-2*Y00[i-1,j]+Y00[i-2,j])/h:end:
if sdx*sdxb>=0 then s1x:=1.0 else s1x:=0.0:end:
vx1:=(Y00[i,j]-Y00[i-1,j])/h+0.5*signum(sdx)*s1x*min(abs(sdx),abs(sdxb)):

if i = N-1 then sdxf:=(Y00[i+1,j]-2*Y00[i+1,j]+Y00[i,j])/h:else sdxf:=(Y00[i+2,j]-2*Y00[i+1,j]+Y00[i,j])/h:end:
if sdx*sdxf>=0 then s1x:=1.0 else s1x:=0.0:end:
vx2:=(Y00[i+1,j]-Y00[i,j])/h-0.5*signum(sdx)*s1x*min(abs(sdx),abs(sdxf)):


if j = 2 then sdyb:=(Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-1])/h:else sdyb:=(Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-2])/h:end:
if sdy*sdyb>=0 then s1y:=1.0 else s1y:=0.0:end:
vy1:=(Y00[i,j]-Y00[i,j-1])/h+0.5*signum(sdy)*s1y*min(abs(sdy),abs(sdyb)):

if j = M-1 then sdyf:=(Y00[i,j+1]-2*Y00[i,j+1]+Y00[i,j])/h:else sdyf:=(Y00[i,j+2]-2*Y00[i,j+1]+Y00[i,j])/h:end:
if sdy*sdyf>=0 then s1y:=1.0 else s1y:=0.0:end:
vy2:=(Y00[i,j+1]-Y00[i,j])/h+0.5*signum(sdy)*s1y*min(abs(sdy),abs(sdyf)):

if v0>=0 then
vx1:=max(vx1,0):vx2:=-min(vx2,0): else
vx1:=-min(vx1,0):vx2:=max(vx2,0): end:

if v0>=0 then
vy1:=max(vy1,0):vy2:=-min(vy2,0): else
vy1:=-min(vy1,0):vy2:=max(vy2,0): end:
nx:=sqrt(max(vx1,vx2)^2+max(vy1,vy2)^2):

F0[i,j]:=nx:
od:od:

end proc:
> 

WENO3:=proc(Y00::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,N::integer,M::integer,Nx::integer,v0::float)
local i::integer,j::integer,h::float[8],nx::float[8],vel::float[8],vx::float[8],vy::float[8],phix::float[8],phiy::float[8],phiave::float[8],jj::integer,phixb::float[8],phixf::float[8],phixb2::float[8],phixf2::float[8],vxb::float[8],phiyb::float[8],phiyf::float[8],vxf::float[8],phiyb2::float[8],phiyf2::float[8],vyb::float[8],vyf::float[8],uf::float[8],ub::float[8],vf::float[8],vb::float[8],tt::float[8],vv0::float[8],sd::float[8],sdf::float[8],sdb::float[8],sdx::float[8],sdy::float[8],sdxb::float[8],s1x::float[8],sdxf::float[8],sdyb::float[8],sdyf::float[8],s1y::float[8],vx1::float[8],vx2::float[8],vy1::float[8],vy2::float[8],w1::float[8],w2::float[8],r1::float[8],r2::float[8],alpha::float[8],beta::float[8],e1::float[8];
e1:=1e-6:

h:=1.0/(N-2):

vv0:=0.1:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
for i from 2 to N-1 do for j from 2 to M-1 do
vx:=0.0:vy:=0.0:phix:=0.0:phiy:=0.0:
phix:=(Y00[i+1,j]-Y00[i-1,j])/2/h:
phiy:=(Y00[i,j+1]-Y00[i,j-1])/2/h:

if i = 2 then sdb:=Y00[i,j]-2*Y00[i-1,j]+Y00[i-1,j]: else sdb:=Y00[i,j]-2*Y00[i-1,j]+Y00[i-2,j]:end:
sd:=Y00[i+1,j]-2*Y00[i,j]+Y00[i-1,j]:
if i = N-1 then sdf:=Y00[i,j]-2*Y00[i+1,j]+Y00[i+1,j]: else sdf:=Y00[i+2,j]-2*Y00[i+1,j]+Y00[i,j]:end:
r1:=(e1+sdb^2)/(e1+sd^2):w1:=1/(1+2*r1^2):r2:=(e1+sdf^2)/(e1+sd^2):w2:=1/(1+2*r2^2):
vx1:=phix-0.5*w1/h*(sd-sdb):
vx2:=phix-0.5*w2/h*(sdf-sd):

if j = 2 then sdb:=Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-1]:else sdb:=Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-2]:end:
sd:=Y00[i,j+1]-2*Y00[i,j]+Y00[i,j-1]:
if j = M-1 then sdf:=Y00[i,j]-2*Y00[i,j+1]+Y00[i,j+1]: else sdf:=Y00[i,j+2]-2*Y00[i,j+1]+Y00[i,j]:end:
r1:=(e1+sdb^2)/(e1+sd^2):w1:=1/(1+2*r1^2):r2:=(e1+sdf^2)/(e1+sd^2):w2:=1/(1+2*r2^2):

vy1:=phiy-0.5*w1/h*(sd-sdb):
vy2:=phiy-0.5*w2/h*(sdf-sd):
alpha:=1.0:beta:=1.0:
if v0>=0 then
vx1:=max(vx1,0):vx2:=-min(vx2,0): else
vx1:=-min(vx1,0):vx2:=max(vx2,0): end:

if v0>=0 then
vy1:=max(vy1,0):vy2:=-min(vy2,0): else
vy1:=-min(vy1,0):vy2:=max(vy2,0): end:
nx:=sqrt(max(vx1,vx2)^2+max(vy1,vy2)^2):
F0[i,j]:=nx:
od:od:

end proc:
> 

UpW:=proc(Y00::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,N::integer,M::integer,Nx::integer,v0::float)
local i::integer,j::integer,h::float[8],nx::float[8],vel::float[8],vx::float[8],vy::float[8],phix::float[8],phiy::float[8],phiave::float[8],jj::integer,phixb::float[8],phixf::float[8],phixb2::float[8],phixf2::float[8],vxb::float[8],phiyb::float[8],phiyf::float[8],vxf::float[8],phiyb2::float[8],phiyf2::float[8],vyb::float[8],vyf::float[8],uf::float[8],ub::float[8],vf::float[8],vb::float[8],tt::float[8],vv0::float[8],sd::float[8],sdf::float[8],sdb::float[8],sdx::float[8],sdy::float[8],sdxb::float[8],s1x::float[8],sdxf::float[8],sdyb::float[8],sdyf::float[8],s1y::float[8],vx1::float[8],vx2::float[8],vy1::float[8],vy2::float[8],e1::float[8],w1::float[8],w2::float[8],r1::float[8],r2::float[8],alpha::float[8],beta::float[8];
e1:=1e-15:

h:=1.0/(N-2):

vv0:=0.1:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
for i from 2 to N-1 do for j from 2 to M-1 do
vx:=0.0:vy:=0.0:phix:=0.0:phiy:=0.0:


vx1:=(Y00[i,j]-Y00[i-1,j])/h:
vx2:=(Y00[i+1,j]-Y00[i,j])/h:


vy1:=(Y00[i,j]-Y00[i,j-1])/h:
vy2:=(Y00[i,j+1]-Y00[i,j])/h:

if v0>=0 then
vx1:=max(vx1,0):vx2:=-min(vx2,0): else
vx1:=-min(vx1,0):vx2:=max(vx2,0): end:

if v0>=0 then
vy1:=max(vy1,0):vy2:=-min(vy2,0): else
vy1:=-min(vy1,0):vy2:=max(vy2,0): end:
nx:=sqrt(max(vx1,vx2)^2+max(vy1,vy2)^2):

F0[i,j]:=nx:
od:od:

end proc:
> 

ENO2:=Compiler:-Compile(ENO2):
> 

WENO3:=Compiler:-Compile(WENO3):
> 

UpW:=Compiler:-Compile(UpW):
> 

printlevel:=1:
> 

 

> 

pFactor := define_external("hw_PardisoFactor", ':-MAPLE', ':-LIB' = "linalg"):
pSolve := define_external("hw_PardisoSolve", ':-MAPLE', ':-LIB' = "linalg"):
pCSF := define_external("hw_SpCompressedSparseForm", ':-MAPLE', ':-LIB' = "linalg"):
> 

#NN:=20;
> 

ss:=proc(NN,delPhi,nn,tf,beta,e1)
local tt,V,V2,V1,TT,dt,vv0,ii,nt,Nt,N,h,w,Nx,MM,M,Ny,ntot,Ntot,Nsp,Y00,Phi0,F0,ff,db,CB,CB0,Ymid,R0,x0,handle,Phiave;
N:=NN+2:h:=1.0/NN:w:=h/2:Nx:=NN:
> 

MM:=NN:M:=MM+2:Ny:=MM:ntot:=N*(M):Ntot:=ntot:Nsp:=5*(N-2)^2+8*(N-1):
> 

Y00:=Matrix(1..N,1..M,datatype=float[8]):F0:=Matrix(1..N,1..M,datatype=float[8]):Phi0:=Vector(1..Ntot,datatype=float[8]):ff:=copy(Phi0):db:=copy(ff):CB0:=Vector(N*M+1,datatype=integer[4]):Ymid:=Matrix(1..N,1..M,datatype=float[8]):
> 

y0proc(NN,Y00,beta,e1):
> 

R0:=Vector(Nsp,datatype=integer[4]):x0:=Vector(Nsp,datatype=float[8]):
> 

CBproc(N,M,CB0):
> 

CB0:=convert(CB0,Vector,datatype=integer[8]):
> 

Rproc(N,M,R0):
> 

R0:=convert(R0,Vector,datatype=integer[8]):
> 

evalf(Eqs11(N,M,Y00,F0,Phi0,evalf(0.1),ff));
> 

XX(N,M,Y00,Phi0,0.1,x0);handle := pFactor(CB0, R0, x0, 11):#handle := pFactor(CB, R, X, 11):
pSolve(ff,db,handle):
> 

PhiAdd(Ntot,Phi0,db):
> 

V2[0]:=(Phi0[ntot-2*N+N/2]/2+Phi0[ntot-2*N+N/2+1]/2)+h/2*0.1;
V1[0]:=Phi0[2*N];
V[0]:=Phi0[ntot-N]+0.5*h*0.1;
> 

TT[0]:=0;tt:=0;
vv0:=max(abs(Phi0[ntot-2*N+1]),abs(Phi0[ntot-2*N+N/2]),abs(Phi0[ntot-2*N+N/2+1]),abs(Phi0[ntot-N])):
> 

dt:=min(h/vv0/nn,tf-tt);
> 

#phiaveAdd(Y00,N,M,Ny);
> 

Phiave[0]:=phiaveAdd(Y00,N,M,Ny);
> 

Nt:=round(tf/dt);
#dt:=tf/Nt;
> 

#Ydat:=Array(1..Nt+1,1..N,1..M,datatype=float[8]):
> 

#YdatStore(Y00,Ydat,N,M,1);
> 

ii:=0:TT[0];tt;
> 

while tt<tf do
delPhi(Y00,Phi0,F0,evalf(dt),N,M,Nx,0.1):
EFAdd(Y00,Ymid,Phi0,F0,dt,N,M,e1);
Eqs11(N,M,Ymid,F0,Phi0,0.1,ff);
XX(N,M,Ymid,Phi0,0.1,x0);handle := pFactor(CB0, R0, x0, 11):
pSolve(ff,db,handle):
PhiAdd(Ntot,Phi0,db):
delPhi(Ymid,Phi0,F0,evalf(dt),N,M,Nx,0.1):
EFAdd2(Y00,Ymid,Phi0,F0,dt,N,M,e1);
Eqs11(N,M,Ymid,F0,Phi0,evalf(0.1),ff);
XX(N,M,Ymid,Phi0,0.1,x0);handle := pFactor(CB0, R0, x0, 11):#handle := pFactor(CB, R, X, 11):
pSolve(ff,db,handle):
PhiAdd(Ntot,Phi0,db):
delPhi(Ymid,Phi0,F0,evalf(dt),N,M,Nx,0.1):
EFAdd3(Y00,Ymid,Phi0,F0,dt,N,M,e1);
Eqs11(N,M,Y00,F0,Phi0,evalf(0.1),ff);
XX(N,M,Y00,Phi0,0.1,x0);handle := pFactor(CB0, R0, x0, 11):#handle := pFactor(CB, R, X, 11):
pSolve(ff,db,handle):
PhiAdd(Ntot,Phi0,db):
ii:=ii+1:#print(i);
V2[ii]:=(Phi0[ntot-2*N+N/2]/2+Phi0[ntot-2*N+N/2+1]/2)+h/2*0.1;#print(ii,V[ii]);
#V1[ii]:=Phi0[ntot]+h/2*0.1;
V1[ii]:=Phi0[2*N];
V[ii]:=Phi0[ntot-N]+0.5*h*0.1;
TT[ii]:=TT[ii-1]+dt:tt:=tt+dt:
#YdatStore(Y00,Ydat,N,M,ii+1);
#Phiave[ii]:=phiaveAdd(Y00,N,M,Ny);
vv0:=max(abs(Phi0[ntot-2*N+1]),abs(Phi0[ntot-2*N+N/2]),abs(Phi0[ntot-2*N+N/2+1]),abs(Phi0[ntot-N])):
#dt:=h/vv0/nn;
dt:=min(h/vv0/nn,tf-tt);#print(tt,ii,dt);
gc();
end:
> 

V2[ii];
> 


end proc:
#time[real]()-t11;time()-t12;NN;
> 

ssapp:=proc(NN,delPhi,nn,tf,beta,e1)
local tt,V,V2,V1,TT,dt,vv0,ii,nt,Nt,N,h,w,Nx,MM,M,Ny,ntot,Ntot,Nsp,Y00,Phi0,F0,ff,db,CB,CB0,Ymid,R0,x0,handle,Phiave, Phimid,Phitemp;
N:=NN+2:h:=1.0/NN:w:=h/2:Nx:=NN:
> 

MM:=NN:M:=MM+2:Ny:=MM:ntot:=N*(M):Ntot:=ntot:Nsp:=5*(N-2)^2+8*(N-1):
> 

Y00:=Matrix(1..N,1..M,datatype=float[8]):F0:=Matrix(1..N,1..M,datatype=float[8]):Phi0:=Vector(1..Ntot,datatype=float[8]):ff:=copy(Phi0):db:=copy(ff):CB0:=Vector(N*M+1,datatype=integer[4]):Ymid:=Matrix(1..N,1..M,datatype=float[8]):Phitemp:=Vector(1..Ntot,datatype=float[8]):Phimid:=Vector(1..Ntot,datatype=float[8]):
> 

y0proc(NN,Y00,beta,e1):
> 

R0:=Vector(Nsp,datatype=integer[4]):x0:=Vector(Nsp,datatype=float[8]):
> 

CBproc(N,M,CB0):
> 

CB0:=convert(CB0,Vector,datatype=integer[8]):
> 

Rproc(N,M,R0):
> 

R0:=convert(R0,Vector,datatype=integer[8]):
> 

evalf(Eqs11(N,M,Y00,F0,Phi0,evalf(0.1),ff));
> 

XX(N,M,Y00,Phi0,0.1,x0);handle := pFactor(CB0, R0, x0, 11):#handle := pFactor(CB, R, X, 11):
pSolve(ff,db,handle):
> 

PhiAdd(Ntot,Phi0,db):
> 

V2[0]:=(Phi0[ntot-2*N+N/2]/2+Phi0[ntot-2*N+N/2+1]/2)+h/2*0.1;
V1[0]:=Phi0[2*N];
V[0]:=Phi0[ntot-N]+0.5*h*0.1;
> 

TT[0]:=0;tt:=0;
vv0:=max(abs(Phi0[ntot-2*N+1]),abs(Phi0[ntot-2*N+N/2]),abs(Phi0[ntot-2*N+N/2+1]),abs(Phi0[ntot-N])):
> 

dt:=h/vv0/nn;
> 

#phiaveAdd(Y00,N,M,Ny);
> 

Phiave[0]:=phiaveAdd(Y00,N,M,Ny);
> 

Nt:=round(tf/dt);
#dt:=tf/Nt;
> 

#Ydat:=Array(1..Nt+1,1..N,1..M,datatype=float[8]):
> 

#YdatStore(Y00,Ydat,N,M,1);
> 

ii:=0:TT[0];tt;
> 

while tt<tf do
PhiSwitch(Ntot,Phi0,Phitemp):
delPhi(Y00,Phi0,F0,evalf(dt),N,M,Nx,0.1):
EFAdd(Y00,Ymid,Phi0,F0,dt,N,M,e1);
Eqs11(N,M,Ymid,F0,Phi0,0.1,ff);
XX(N,M,Ymid,Phi0,0.1,x0);handle := pFactor(CB0, R0, x0, 11):
pSolve(ff,db,handle):
PhiAdd(Ntot,Phi0,db):
MidPred(Ntot,Phi0,Phitemp,Phimid):
delPhi(Ymid,Phi0,F0,evalf(dt),N,M,Nx,0.1):
EFAdd2(Y00,Ymid,Phi0,F0,dt,N,M,e1);
delPhi(Ymid,Phimid,F0,evalf(dt),N,M,Nx,0.1):
EFAdd3(Y00,Ymid,Phimid,F0,dt,N,M,e1);
ii:=ii+1:#print(i);
V2[ii]:=(Phi0[ntot-2*N+N/2]/2+Phi0[ntot-2*N+N/2+1]/2)+h/2*0.1;#print(ii,V[ii]);
#V1[ii]:=Phi0[ntot]+h/2*0.1;
V1[ii]:=Phi0[2*N];
V[ii]:=Phi0[ntot-N]+0.5*h*0.1;
TT[ii]:=TT[ii-1]+dt:tt:=tt+dt:
#YdatStore(Y00,Ydat,N,M,ii+1);
#Phiave[ii]:=phiaveAdd(Y00,N,M,Ny);
vv0:=max(abs(Phi0[ntot-2*N+1]),abs(Phi0[ntot-2*N+N/2]),abs(Phi0[ntot-2*N+N/2+1]),abs(Phi0[ntot-N])):
#dt:=h/vv0/nn;
dt:=min(h/vv0,tf-tt);#print(tt,ii,dt);
gc();
end:
V2[ii];
end proc:
#time[real]()-t11;time()-t12;NN;
> 

 

A procedure is written to find the value of the potential at t =tf using different upwind methods usinb both SSP RK3 and SSP RK3a. Nt1 can be set to 8 to get results till 1280.

> 

Nt1:=5;

5

 (4)
> 

for i from 1 to Nt1 do
ntot:=10*2^(i-1):
t11:=time[real]():t12:=time():ss1[i]:=ss(ntot,WENO3,1,2.0,0.5,1e-9):
if i=1 then err:=1: else err:=ss1[i]-ss1[i-1]:end:
print(time()-t12,time[real]()-t11,ntot,err,ss1[i]):od:

3.954, .687, 10, 1, HFloat(0.1682134687160814)

  

6.359, .814, 20, HFloat(0.0010155538425478672), HFloat(0.16922902255862926)

  

14.922, 1.911, 40, HFloat(6.930403646316641e-4), HFloat(0.16992206292326092)

  

38.250, 4.858, 80, HFloat(2.4869671206081967e-4), HFloat(0.17017075963532174)

  

169.265, 21.364, 160, HFloat(8.121914997127888e-5), HFloat(0.17025197878529302)

 (5)
> 

for i from 1 to Nt1 do
ntot:=10*2^(i-1):
t11:=time[real]():t12:=time():ss1[i]:=ss(ntot,ENO2,1,2.0,0.5,1e-9):
if i=1 then err:=1: else err:=ss1[i]-ss1[i-1]:end:
print(time()-t12,time[real]()-t11,ntot,err,ss1[i]):od:

3.015, .378, 10, 1, HFloat(0.1671674213285358)

  

6.547, .827, 20, HFloat(0.0015764859305227308), HFloat(0.16874390725905852)

  

15.094, 1.947, 40, HFloat(9.851888005965037e-4), HFloat(0.16972909605965503)

  

38.547, 4.925, 80, HFloat(3.6093200578019013e-4), HFloat(0.17009002806543522)

  

171.578, 21.645, 160, HFloat(1.251913120496606e-4), HFloat(0.17021521937748488)

 (6)
> 

for i from 1 to Nt1 do
ntot:=10*2^(i-1):
t11:=time[real]():t12:=time():ss1[i]:=ss(ntot,UpW,1,2.0,0.5,1e-9):
if i=1 then err:=1: else err:=ss1[i]-ss1[i-1]:end:
print(time()-t12,time[real]()-t11,ntot,err,ss1[i]):od:

3.000, .375, 10, 1, HFloat(0.16652718190712018)

  

6.375, .807, 20, HFloat(5.646156758065679e-5), HFloat(0.16658364347470084)

  

15.000, 1.935, 40, HFloat(0.0015059794761171363), HFloat(0.16808962295081797)

  

37.703, 4.811, 80, HFloat(0.001027021876081169), HFloat(0.16911664482689914)

  

163.578, 20.624, 160, HFloat(5.517479717490614e-4), HFloat(0.1696683927986482)

 (7)
> 

Nt1:=5;

5

 (8)
> 

for i from 1 to Nt1 do
ntot:=10*2^(i-1):
t11:=time[real]():t12:=time():ss1[i]:=ssapp(ntot,WENO3,1,2.0,0.5,1e-9):
if i=1 then err:=1: else err:=ss1[i]-ss1[i-1]:end:
print(time()-t12,time[real]()-t11,ntot,err,ss1[i]):od:

3.375, .416, 10, 1, HFloat(0.16807139153530765)

  

6.329, .809, 20, HFloat(0.0011618505716259886), HFloat(0.16923324210693363)

  

12.218, 1.567, 40, HFloat(7.037260277996393e-4), HFloat(0.16993696813473327)

  

29.672, 3.805, 80, HFloat(2.4502719840738263e-4), HFloat(0.17018199533314066)

  

90.485, 11.507, 160, HFloat(7.722404563609286e-5), HFloat(0.17025921937877675)

 (9)
> 

 

> 

for i from 1 to Nt1 do
ntot:=10*2^(i-1):
t11:=time[real]():t12:=time():ss1[i]:=ssapp(ntot,ENO2,1,2.0,0.5,1e-9):
if i=1 then err:=1: else err:=ss1[i]-ss1[i-1]:end:
print(time()-t12,time[real]()-t11,ntot,err,ss1[i]):od:

3.266, .422, 10, 1, HFloat(0.16713029032349638)

  

6.500, .832, 20, HFloat(0.0016324550833170803), HFloat(0.16876274540681346)

  

12.578, 1.603, 40, HFloat(9.833398212751587e-4), HFloat(0.16974608522808862)

  

28.656, 3.662, 80, HFloat(3.5456335900937996e-4), HFloat(0.170100648587098)

  

95.203, 12.123, 160, HFloat(1.2050296782994963e-4), HFloat(0.17022115155492795)

 (10)
> 

 

> 

 

> 

for i from 1 to Nt1 do
ntot:=10*2^(i-1):
t11:=time[real]():t12:=time():ss1[i]:=ssapp(ntot,UpW,1,2.0,0.5,1e-9):
if i=1 then err:=1: else err:=ss1[i]-ss1[i-1]:end:
print(time()-t12,time[real]()-t11,ntot,err,ss1[i]):od:

3.125, .407, 10, 1, HFloat(0.16629594073351656)

  

6.343, .820, 20, HFloat(2.919030713841586e-4), HFloat(0.16658784380490071)

  

12.875, 1.635, 40, HFloat(0.0015038491873215765), HFloat(0.1680916929922223)

  

29.375, 3.770, 80, HFloat(0.0010269012198876326), HFloat(0.16911859421210992)

  

87.282, 11.049, 160, HFloat(5.512458439223822e-4), HFloat(0.1696698400560323)

 (11)
> 

 


 

Download Pardisocodetoupload.mw

 

any takers?...

Dr. Venkat Subramanian 386
March 07 2022
0 0

Just checking to see if this can be done. We are in the process of writing a paper, and this can help in speeding up the code shared with the community.

Great question...

Dr. Venkat Subramanian 386
February 07 2022
0 0

Unfortunately, you may not get much help on this or general questions about DLL (for Windows). 

Instead of using many free solvers (ODES/DAEs) in C/Fortran and calling them as DLL, Maple, unfortunately, moved to rewrite inbuilt versions which work well for many situations, but will always lag behind the C codes.

It may be the compiler issue, but in general, I haven't been successful in easily integrating multiple C codes in DLL form. You may be better off writing to C from Maple and then executing the C file from Maple using a system command if you want efficiency.

PS, it is possible to write procedures in Maple and compile the same for efficiency. But storing compiled codes and reusing them is not trivial and involves the use of undocumented commands.

numerical solution is attached...

Dr. Venkat Subramanian 386
September 19 2021
0 0

z goes only from - 1 to 0. 

u1,theta1,N solved from -1 to 0.
u2, theta2 are also solved from -1 to 0 (for z)
I think you have 3 variables in region 1 and 2 in region 2. variables and fluxes are continuous at y = 0. 

it will be great in the future if you state all the bcs at y = -1 first, then at y =0, and then at y = 1, for the convenience of coding. 

Until Maple 17 or so D(u[1])(0) worked. In the later version, you have to use u1.
 

> 

restart;

BVP is solved in z. z varies from -1 to 0.

> 

eq1:=diff(u1(z),z,z) + 1/2*diff(N(z),z) + 1/2*theta1(z) = 0;
 

eq1 := diff(u1(z), `$`(z, 2))+(1/2)*(diff(N(z), z))+(1/2)*theta1(z) = 0

 (1)
> 

eq2:=diff(N(z),z$2) - 2/3 * (2*N(z) + diff(u1(z),z)) = 0;
 

eq2 := diff(N(z), `$`(z, 2))-(4/3)*N(z)-(2/3)*(diff(u1(z), z)) = 0

 (2)
> 

eq3:=diff(theta1(z),z$2) = 0;

eq3 := diff(theta1(z), `$`(z, 2)) = 0

 (3)
> 

eq4:=diff(u2(z),z$2) + theta2(z) = 0;

eq4 := diff(u2(z), `$`(z, 2))+theta2(z) = 0

 (4)
> 

eq5:=diff(theta2(z),z$2) = 0;

eq5 := diff(theta2(z), `$`(z, 2)) = 0

 (5)
> 

bcs:=u1(-1)=0,theta1(-1)=1,N(-1)=0,u2(0)=0,theta2(0)=0,u1(0)=u2(-1),D(N)(0)=0,theta1(0)=theta2(-1),D(theta1)(0)=D(theta2)(-1),D(u1)(0)+N(0)/2=1/2*D(u2)(-1);

bcs := u1(-1) = 0, theta1(-1) = 1, N(-1) = 0, u2(0) = 0, theta2(0) = 0, u1(0) = u2(-1), (D(N))(0) = 0, theta1(0) = theta2(-1), (D(theta1))(0) = (D(theta2))(-1), (D(u1))(0)+(1/2)*N(0) = (1/2)*(D(u2))(-1)

 (6)
> 

soln:=dsolve({eq1,eq2,eq3,eq4,eq5,bcs},type=numeric):
> 

soln(-1);

[z = -1., N(z) = 0., diff(N(z), z) = -0.904317007826591e-1, theta1(z) = 1.00000000000000, diff(theta1(z), z) = -.500000000000000, theta2(z) = .500000000000000, diff(theta2(z), z) = -.500000000000000, u1(z) = 0., diff(u1(z), z) = .371898088450732, u2(z) = .172870489765202, diff(u2(z), z) = -0.620382309853562e-2]

 (7)
> 

soln(0);

[z = 0., N(z) = -0.228483654680803e-1, diff(N(z), z) = 0., theta1(z) = .500000000000000, diff(theta1(z), z) = -.500000000000000, theta2(z) = 0., diff(theta2(z), z) = -.500000000000000, u1(z) = .172870489765202, diff(u1(z), z) = 0.832227118477233e-2, u2(z) = 0., diff(u2(z), z) = -.256203823098536]

 (8)
> 

p1:=plots:-odeplot(soln,[z,u1(z)],-1..0):
> 

p2:=plots:-odeplot(soln,[z+1,u2(z)],-1..0,color=green):
> 

plots:-display({p1,p2},axes=boxed,labels=[y,"u1,u2"]);

 
> 

 


Please confirm that this is ok.

Download threepointBVPmodify.mws

Why not vectors?...

Dr. Venkat Subramanian 386
July 19 2021
0 0

@acer 

Since this example has only two inputs, f converts the vector to scalars.

It will be great to have fdiff which takes obj which takes a vector and performs the calculation. When we have 10 or 20 inputs to the objective, it is a pain to convert the model to scalars from vectors. In that case, why can't we write obj in the operator form and use it directly?

 

stiff =true...

Dr. Venkat Subramanian 386
April 30 2021
0 0

@tomleslie 

Please note that stiff error comes when stiff=true is used. This is the failure of Rosenbrock methods, probably from Maple's failure to find the analytic Jacobian.

For the OP, you may be better off writing your system as a system of first order equations (Maple does this before going to dsolve). Also try implicit =true to see if that helps. This can be done only for first order ODEs.

NLPSolve should be modified to take care...

Dr. Venkat Subramanian 386
March 23 2021
0 0

@tomleslie 

Thanks for your solution. I believe the NLPSolve with method=nonlinearsimplex has no issues (it can't deal with bounds and constraints). It would be great if NLPSolve can use an inbuilt wrapper to make the optimization work for simple returns from the procedures without these _passed numeric returns. I believe the issue may be due to the calls for fdiff. 
Please note that GlobalOptimzation package calls the sqp optimizer, and has no issues with procedures without the extra _passed commands. 

update...

Dr. Venkat Subramanian 386
March 14 2021

@ I am able to run classicworksheet with maple2021

I had to go to to bin.win folder

run cwmaple.exe

This is actually better than previous versions where I could only run 32bit version in classic worksheet mode. The help file says that cwmaple.exe exists inside bin.x86 folder,  but it is in bin.win folder

till 2019...

Dr. Venkat Subramanian 386
March 11 2021

@nm 

I have classic worksheet working for Maple2019 (a little unstable, but it works). So classic worksheet (32bit) was dropped either in 2020 or 2021.

Yet to check the additions...

Dr. Venkat Subramanian 386
March 11 2021

Unfortunately no more classic worksheet. I understand that there may not be a demand for 32 bit installations, but I am a big fan of classic worksheet interface, please add it if possible.

Not true...

Dr. Venkat Subramanian 386
January 20 2021
0 0

@tomleslie 

Please note that this statement (IVP numerical methods are more stable compared to BVP methods) is incorrect. Actually, BVP methods are more stable to solve a given BVP problem. It is easy to show many examples to prove the same. For example, 

take dy1/dt = y2

dy2/dt = tau^2*y1 + (Pi^2-tau^2)*sin(Pi*t)

IVP: t=0, y1= 0, y2 = Pi

BVP: t= 0, y1 = 0; t = 1, y1 = 0

(Credit: Dr. Biegler at CMU)

This model is very stable with a BVP approach, but not with IVP approach. Shooting method is not useful for largescale or stiff problems. At least multiple shooting is needed.

Add some noise to ics, constants, you will see BVP methods behaving much better.

 

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