Venkat Subramanian

386 Reputation

13 Badges

14 years, 165 days

MaplePrimes Activity


These are questions asked by

@acer,
In the past, you had shown how to use take a given sparse A matrix and store the LU components and use the same factorization for multiple backsolves at https://www.mapleprimes.com/posts/41191-Solving-Sparse-Linear-Systems-In-Maple#comment200817

When PDEs are solved, we are calling the A matrix at every time step (with some discretizations in x) and the factorization is done and stored at every time step (this is time-consuming and memory-consuming). Is it possible to call UMFPACK using only the sparse storage and entries? The main routine seems to provide the option 

https://people.sc.fsu.edu/~jburkardt/f77_src/umfpack/umfpack.html

This would mean that the pattern is found only once for the Jacobian at t= 0, then only the non-zero sparse entries (vector/row, not a matrix) are updated at every time step.

To be clear, what I am asking for is create a random sparse matrix (say 4x4 or 10x10)
create a b Vector.

Solve AZ =b with method = SparseDirect. (This should internally create and store R, CC, X as at https://people.sc.fsu.edu/~jburkardt/f77_src/umfpack/umfpack.html)

Store R, CC, X and update only X for a different A matrix with the same sparsity pattern and solve for Z using stored (R, CC), and updated X.

Hi All,

https://www.maplesoft.com/support/helpJP/Maple/view.aspx?path=NAG%2ff11dec

This link shows that there are 3 options rgmres, cgs and bicgstab. When I use iterative solver, it uses the cgs option.

LinearAlgebra:-LinearSolve(
      A, b, method= SparseIterative)

LinearSolve:   "copying right hand side to enable external call"
LinearSolve:   "using method"   SparseIterative
LinearSolve:   "using method  "   SparseIterative
LinearSolve:   "calling external function"
LinearSolve:   "using CGS method"
LinearSolve:   "preconditioning with incomplete LU factorization"
LinearSolve:   "level of fill = "   0
LinearSolve:   "using complete pivoting strategy"
LinearSolve:   "dimension of workspaces for preconditioner = "   44
LinearSolve:   "using infinity norm in stopping criteria"
LinearSolve:   "setting maximum iterations to "   200
LinearSolve:   "setting tolerance to "   .10e-7
LinearSolve:   "NAG"   hw_f11zaf
LinearSolve:   "NAG"   hw_f11daf
LinearSolve:   "NAG"   hw_f11dcf
LinearSolve:   "number of iterations"   1
LinearSolve:   "residual computed last as"   3.33066907387547e-016
 

 

How can I force Maple to use BiCGSTAB?

thanks

I understand how Browse() works. Sometimes Maple crashes or fails for some programs. Is it possible to add print("here" ) into the the mla files in the library to identify why and where Maple fails?

 

Thanks

 

I would like to return local variable y (line 4 in showstat) in the attached dummy procedure (s1) without manually adding any comment inside the procedure s1. This procedure is a simple one and easy to copy paste/or change. When we have a long procedure, it is difficult to do so. I will always know the name of the local variable I want (say, y) and/or line number in showstat

Thanks

PS, I want to get y:=array(1..,2[(1)=x,2=zz])
 

restart;

s1:=proc(n,x)
local y,xx,i,j,zz::array(1..n,1..n);
for i from 1 to n do for j from 1 to n do zz[i,j]:=x[i]*(1+x[j]^2);od:
od:
y:=array(1..2,[(1)=x, (2)=zz]):
for j from 1 to n do xx[i]:=zz[i,i]/(add(zz[i,j],j=1..n));od:
0;
end proc;

s1 := proc (n, x) local y, xx, i, j, zz::(array(1 .. n, 1 .. n)); for i to n do for j to n do zz[i, j] := x[i]*(1+x[j]^2) end do end do; y := array(1 .. 2, [1 = x, 2 = zz]); for j to n do xx[i] := zz[i, i]/add(zz[i, j], j = 1 .. n) end do; 0 end proc

(1)

showstat(s1);

 

s1 := proc(n, x)

local y, xx, i, j, zz::array(1 .. n,1 .. n);

   1   for i to n do

   2     for j to n do

   3       zz[i,j] := x[i]*(1+x[j]^2)

         end do

       end do;

   4   y := array(1 .. 2,[1 = x, 2 = zz]);

   5   for j to n do

   6     xx[i] := zz[i,i]/add(zz[i,j],j = 1 .. n)

       end do;

   7   0

end proc

 

 

x0:=Vector(2,[1,1]);

x0 := Vector(2, {(1) = 1, (2) = 1})

(2)

s1(2,x0);

0

(3)

 


 

Download showstatexample.mws


To summarize,

dsolve numeric parametric form gives wrong answers when fdiff is used to calculate the Hessian. Doing dsolve twice gives the correct answer. Hope I am not making syntax or programming errors.

Also, fdiff is not compatible with vector form.

restart;

Digits:=15;

Digits := 15

(1)

sys:=diff(y1(t),t)=-(u+u^2/2)*y1(t),diff(y2(t),t)=u*y1(t);

sys := diff(y1(t), t) = -(u+(1/2)*(u^2))*y1(t), diff(y2(t), t) = u*y1(t)

(2)

 The system is solved with initial conditions 1,0 and value of u being equal to u1 for t<=0.5, and u2 for t>0.5. The objective is y2 at t =1. There are errors when Hessian is calculated using parametric dsolve and fdiff

sol1 := dsolve({sys, y1(0) = alpha, y2(0) = beta}, type = numeric, 'parameters' = [alpha, beta, u],maxfun=0,range=0..0.5):

obj1:=proc(u1,u2)
local z1,s1,s2,z2;
global sol1;
sol1('parameters'=[1,0,u1]);
z1 := sol1(.5);
s1:=subs(z1,y1(t));
s2:=subs(z1,y2(t));
sol1('parameters'=[s1,s2,u2]):
z2:=sol1(0.5);
-subs(z2,y2(t));
end proc;

obj1 := proc (u1, u2) local z1, s1, s2, z2; global sol1; sol1('parameters' = [1, 0, u1]); z1 := sol1(.5); s1 := subs(z1, y1(t)); s2 := subs(z1, y2(t)); sol1('parameters' = [s1, s2, u2]); z2 := sol1(.5); -subs(z2, y2(t)) end proc

(3)

u0:=[0.8,1.8];

u0 := [.8, 1.8]

(4)

obj1(op(u0));

-.552540796143903

(5)

Hess1:=Matrix(2,2):for i from 1 to 2 do for j from 1 to 2 do Hess1[i,j]:=fdiff(obj1,[i,j],u0,workprec=1.0);od;od;Hess1;

Matrix(2, 2, {(1, 1) = 115.470000000000, (1, 2) = -1.11000000000000, (2, 1) = -1.11000000000000, (2, 2) = 2.22000000000000})

(6)

 Hessian is wrong. It should be noted that when dsolve is called twice inside the objective function there is no error in Hessian calculation with fdiff. The parametric dsolve is efficient, how to get correct answers with that? Another question is how to use fdiff with vector form instead u1,u2 etc to make it generic for any number of variables?

obj2:=proc(u1,u2)
local z1,s1,s2,z2,sol2;
sol2:=dsolve({op(subs(u=u1,[sys])), y1(0) = 1, y2(0) = 0}, type = numeric,maxfun=0,range=0..0.5);
z1 := sol2(.5);
s1:=subs(z1,y1(t));
s2:=subs(z1,y2(t));
sol2:=dsolve({op(subs(u=u2,[sys])), y1(0) = s1, y2(0) = s2}, type = numeric,maxfun=0,range=0..0.5):
z2:=sol2(0.5);
-subs(z2,y2(t));
end proc;

obj2 := proc (u1, u2) local z1, s1, s2, z2, sol2; sol2 := dsolve({op(subs(u = u1, [sys])), y1(0) = 1, y2(0) = 0}, type = numeric, maxfun = 0, range = 0 .. .5); z1 := sol2(.5); s1 := subs(z1, y1(t)); s2 := subs(z1, y2(t)); sol2 := dsolve({op(subs(u = u2, [sys])), y1(0) = s1, y2(0) = s2}, type = numeric, maxfun = 0, range = 0 .. .5); z2 := sol2(.5); -subs(z2, y2(t)) end proc

(7)

obj2(op(u0));

-.552540796143903

(8)

Hess2:=Matrix(2,2):for i from 1 to 2 do for j from 1 to 2 do Hess2[i,j]:=fdiff(obj2,[i,j],u0,workprec=1.0);od;od;Hess2;

Matrix(2, 2, {(1, 1) = .234295697631809, (1, 2) = 0.101868715004724e-1, (2, 1) = 0.101868715004724e-1, (2, 2) = 0.853480932660994e-1})

(9)

 Doing dsolve twice gives correct answer compared to one parametric dsolve even though objective returns the same number (even gradient returns the expected answers).


 

Download Bugindsolveorfdiff.mws

1 2 3 Page 2 of 3