## 406 Reputation

14 years, 320 days

## Bug in dsolve numeric or fdiff?...

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;
 (1)
 > sys:=diff(y1(t),t)=-(u+u^2/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;
 (3)
 > u0:=[0.8,1.8];
 (4)
 > obj1(op(u0));
 (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;
 (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;
 (7)
 > obj2(op(u0));
 (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;
 (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).

## Maple 6 version...

Maple

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).

## Towards a compiled Newton Raphson...

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;
 (1)
 > N:=4;
 (2)
 > f:=array([seq(evalf(-2*x[i]+x[i]^2+0.99/i^2),i=1..N)]);
 (3)
 > fsolve({f[1],f[2],f[3],f[4]},{x[1]=1.2,x[2]=1,x[3]=0,x[4]=0});
 (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]);
 (5)
 > f0:=Vector([seq(0,i=1..N)],datatype=float[8]);
 (6)
 > eqsJ := [seq(seq(Jac[i,j]=evalf(diff(f[i],x[j])),j=1..N),i=1..N)];
 (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]);
 (8)
 > f0;
 (9)
 > j0:=Matrix(1..N,1..N,datatype=float[8]);
 (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;
 >
 (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;
 (13)