Venkat Subramanian

386 Reputation

13 Badges

14 years, 168 days

MaplePrimes Activity


These are questions asked by

I have been in touch with Maplesoft trying to get this version for windows (they are not able to create a download for this). I have codes that used to run in Maple 6 but not in Maple 7 or later. (Maple V should work as well).

If you any of have this version, please let me know if I can try it out for a limited time (I have always had licenses from Maple V Release 3 or 4).

I am not able to post those codes for obvious confidentiality reasons.

 

Thanks

 

(I tried my code in Maple 7, but no use). 

Acer, 

Thanks. This request from me comes from the need that Maple does not have a real DAE solver. (It converts ODEs to DAEs and then solves, which fails for nonlinear and most non-trivial DAEs).

Your question is valid. For situations which can be compiled (i.e, no Heaviside or functions that can't be compiled) at hardware floats, I would like to have a Newton Raphson solver. Some parameters (time values) change when you integrate DAEs, which means that there is a need to solve the same set of nonlinear equations again and again for minor changes in parameters.

While LinearSolve is fast, evaluating F and Jac for a large system of equations can kill the speed. I am attaching an example for Newton wherein F, Jac and a linear solver is compiled. But the Newton procedure is not. If we can compile the Newton procedure also, it will be much more faster.

Also, having a compiled code with parametric form will help in optimization. 

restart:

Digits:=15;

Digits := 15

(1)

N:=4;

N := 4

(2)

f:=array([seq(evalf(-2*x[i]+x[i]^2+0.99/i^2),i=1..N)]);

f := Vector[row](4, {(1) = -2.*x[1]+x[1]^2+.99, (2) = -2.*x[2]+x[2]^2+.247500000000000, (3) = -2.*x[3]+x[3]^2+.110000000000000, (4) = -2.*x[4]+x[4]^2+0.618750000000000e-1})

(3)

fsolve({f[1],f[2],f[3],f[4]},{x[1]=1.2,x[2]=1,x[3]=0,x[4]=0});

{x[1] = .900000000000000, x[2] = .132532421355126, x[3] = 1.94339811320566, x[4] = 0.314314686094742e-1}

(4)

 

eqsA := [seq(F[i]=f[i],i=1..N)]:

irform2 := StatSeq(seq(Assign(op(i)),i=eqsA)):

prccons:= codegen[intrep2maple](Proc(Parameters(x::Array,F::Array),irform2)):

f3:=Compiler:-Compile(prccons):

 

 

x0:=Vector([seq(1/2,i=1..N)],datatype=float[8]);

x0 := Vector(4, {(1) = .500000000000000, (2) = .500000000000000, (3) = .500000000000000, (4) = .500000000000000})

(5)

f0:=Vector([seq(0,i=1..N)],datatype=float[8]);

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

(6)

eqsJ := [seq(seq(Jac[i,j]=evalf(diff(f[i],x[j])),j=1..N),i=1..N)];

eqsJ := [Jac[1, 1] = -2.+2.*x[1], Jac[1, 2] = 0., Jac[1, 3] = 0., Jac[1, 4] = 0., Jac[2, 1] = 0., Jac[2, 2] = -2.+2.*x[2], Jac[2, 3] = 0., Jac[2, 4] = 0., Jac[3, 1] = 0., Jac[3, 2] = 0., Jac[3, 3] = -2.+2.*x[3], Jac[3, 4] = 0., Jac[4, 1] = 0., Jac[4, 2] = 0., Jac[4, 3] = 0., Jac[4, 4] = -2.+2.*x[4]]

(7)

 

irformJ := StatSeq(seq(Assign(op(i)),i=eqsJ)):

prcconsJ:= codegen[intrep2maple](Proc(Parameters(x::Array,Jac::Array),irformJ)):

J3:=Compiler:-Compile(prcconsJ):

 

 

j0:=Matrix(1..N,1..N,datatype=float[8]):

 Following linear sovler algorithm is basiclly from dsolve/numeric/SC

.

 

 

s3:=proc(n::posint, A::Array( datatype = float[ 8 ] ) , ip::Array( datatype = integer[ 4 ] ),b::Array( datatype = float[ 8 ] ) )
local i::integer, j::integer, k::integer, m::integer, t::float;
      ip[n] := 1;
      for k to n-1 do
        m := k;
        for i from k+1 to n do
          if abs(A[m,k]) < abs(A[i,k]) then
            m := i
           end if
         end do;
        ip[k] := m;
        if m <> k then
          ip[n] := -ip[n]
         end if;
       t := A[m,k];
       A[m,k] := A[k,k];
       A[k,k] := t;
       if t = 0 then
         ip[n] := 0;
         #return ip
         end if;
       for i from k+1 to n do
         A[i,k] := -A[i,k]/t
         end do;
       for j from k+1 to n do
         t := A[m,j];
         A[m,j] := A[k,j];
         A[k,j] := t;
         if t <> 0 then
           for i from k+1 to n do
             A[i,j] := A[i,j]+A[i,k]*t
             end do
           end if
         end do
       end do;
     if A[n,n] = 0 then
       ip[n] := 0
       end if;
     #ip[n]
      if ip[n] = 0 then
        error "The matrix A is numerically singular"
       end if;
      if 1 < n then
        for k to n-1 do
          m := ip[k];
          t := b[m];
          b[m] := b[k];
          b[k] := t;
          for i from k+1 to n do
           b[i] := b[i]+A[i,k]*t
           end do
         end do;
       for k from n by -1 to 2 do
         b[k] := b[k]/A[k,k];
         t := -b[k];
        for i to k-1 do
           b[i] := b[i]+A[i,k]*t
           end do
         end do
       end if;
     b[1] := b[1]/A[1,1];
return;
end proc:

s3c:=Compiler:-Compile(s3):

 

 

x0:=Vector(N,[seq(0.5,i=1..N)],datatype=float[8]);

x0 := Vector(4, {(1) = .500000000000000, (2) = .500000000000000, (3) = .500000000000000, (4) = .500000000000000})

(8)

f0;

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

(9)

j0:=Matrix(1..N,1..N,datatype=float[8]);

j0 := Matrix(4, 4, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (2, 1) = 0., (2, 2) = 0., (2, 3) = 0., (2, 4) = 0., (3, 1) = 0., (3, 2) = 0., (3, 3) = 0., (3, 4) = 0., (4, 1) = 0., (4, 2) = 0., (4, 3) = 0., (4, 4) = 0.})

(10)

10

(11)

Newton:=proc(x0,f0,j0,N)
global f3,s3c;
local xnew,err,xold,ip,i::integer;
ip:=Vector(N,[seq(1,i=1..N)],datatype=integer[4]);
err:=10;
xold:=x0;
while err>1e-8 do
f3(xold,f0);
J3(xold,j0);
s3c(N,-j0,ip,f0);
xnew:=xold+f0;
err:=max(abs(f0))/N;
xold:=xnew;
end:
xnew;
end proc;

 

Newton := proc (x0, f0, j0, N) local xnew, err, xold, ip, i::integer; global f3, s3c; ip := Vector(N, [seq(1, i = 1 .. N)], datatype = integer[4]); err := 10; xold := x0; while 0.1e-7 < err do f3(xold, f0); J3(xold, j0); s3c(N, -j0, ip, f0); xnew := xold+f0; err := max(abs(f0))/N; xold := xnew end do; xnew end proc

(12)

Compiler:-Compile(Newton);

Error, (in CodeGeneration:-IssueError) cannot analyze non-integer range boundary args[4]

 

t11:=time():Newton(x0,f0,j0,N);time()-t11;

Vector(4, {(1) = .900000000000000, (2) = .132532421355126, (3) = 0.566018867943396e-1, (4) = 0.314314686094742e-1})

0.

(13)


Download NewtonRaphsoncompileattempt.mws

1 2 3 Page 3 of 3