Advanced Search

Question: Bug in dsolve numeric or fdiff?

Question:Bug in dsolve numeric or fdiff?

Dr. Venkat Subramanian 386
dsolve numeric differential-equations
November 03 2016
0


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
Please Wait...