## Copyright (C) 2008 Carlo de Falco
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, see .
## -*- texinfo -*-
## @deftypefn {Function File} {@var{A}} = bvp4c (@var{odefun}, @var{bcfun}, @var{solinit})
##
## Solves the first order system of non-linear differential equations defined by
## @var{odefun} with the boundary conditions defined in @var{bcfun}.
##
## The structure @var{solinit} defines the grid on which to compute the
## solution (@var{solinit.x}), and an initial guess for the solution (@var{solinit.y}).
## The output @var{sol} is also a structure with the following fields:
## @itemize
## @item @var{sol.x} list of points where the solution is evaluated
## @item @var{sol.y} solution evaluated at the points @var{sol.x}
## @item @var{sol.yp} derivative of the solution evaluated at the
## points @var{sol.x}
## @item @var{sol.solver} = "bvp4c" for compatibility
## @end itemize
## @seealso{odpkg}
## @end deftypefn
## Author: Carlo de Falco
## Created: 2008-09-05
function sol = bvp4cold(odefun,bcfun,solinit,options)
if (isfield(solinit,"x"))
t = solinit.x;
else
error("bvp4c: missing initial mesh solinit.x");
end
if (isfield(solinit,"y"))
u_0 = solinit.y;
else
error("bvp4c: missing initial guess");
end
if (isfield(solinit,"parameters"))
error("bvp4c: solving for unknown parameters is not yet supported");
end
RelTol = 1e-3;
AbsTol = 1e-6;
if ( nargin > 3 )
if (isfield(options,"RelTol"))
RelTol = options.RelTol;
endif
if (isfield(options,"RelTol"))
AbsTol = options.AbsTol;
endif
endif
Nvar = rows(u_0);
Nint = length(t)-1;
s = 3;
h = diff(t);
AbsErr = inf;
RelErr = inf;
MaxIt = 10;
for iter = 1:MaxIt
x = [ u_0(:); zeros(Nvar*Nint*s,1) ];
x = __bvp4c_solve__ (t, x, h, odefun, bcfun, Nvar, Nint, s);
u = reshape(x(1:Nvar*(Nint+1)),Nvar,Nint+1);
for kk=1:Nint+1
du(:,kk) = odefun(t(kk), u(:,kk));
end
tm = (t(1:end-1)+t(2:end))/2;
um = [];
for nn=1:Nvar
um(nn,:) = interp1(t,u(nn,:),tm);
endfor
f_est = [];
for kk=1:Nint
f_est(:,kk) = odefun(tm(kk), um(:,kk));
end
du_est = [];
for nn=1:Nvar
du_est(nn,:) = diff(u(nn,:))./h;
end
err = max(abs(f_est-du_est));
AbsErr = max(err)
RelErr = AbsErr/norm(du,inf)
if ( (AbsErr >= AbsTol) && (RelErr >= RelTol) )
ref_int = find( (err >= AbsTol) & (err./max(max(abs(du))) >= RelTol) );
t_add = tm(ref_int);
t_old = t;
t = sort([t, t_add]);
h = diff(t);
u_0 = [];
for nn=1:Nvar
u_0(nn,:) = interp1(t_old, u(nn,:), t);
end
Nvar = rows(u_0);
Nint = length(t)-1
else
break
end
endfor
## K = reshape(x([1:Nvar*Nint*s]+Nvar*(Nint+1)),Nvar,Nint,s);
## K1 = reshape(K(:,:,1), Nvar, Nint);
## K2 = reshape(K(:,:,2), Nvar, Nint);
## K3 = reshape(K(:,:,3), Nvar, Nint);
sol.x = t;
sol.y = u;
sol.yp= du;
sol.parameters = [];
sol.solver = 'bvp4c';
endfunction
function diff_K = __bvp4c_fun_K__ (t, u, Kin, f, h, s, Nint, Nvar)
%% coefficients
persistent C = [0 1/2 1 ];
persistent A = [0 0 0;
5/24 1/3 -1/24;
1/6 2/3 1/6];
for jj = 1:s
for kk = 1:Nint
Y = repmat(u(:,kk),1,s) + ...
(reshape(Kin(:,kk,:),Nvar,s) * A.') * h(kk);
diff_K(:,kk,jj) = Kin(:,kk,jj) - f (t(kk)+C(jj)*h(kk), Y);
endfor
endfor
endfunction
function diff_u = __bvp4c_fun_u__ (t, u, K, h, s, Nint, Nvar)
%% coefficients
persistent B= [1/6 2/3 1/6 ];
Y = zeros(Nvar, Nint);
for jj = 1:s
Y += B(jj) * K(:,:,jj);
endfor
diff_u = u(:,2:end) - u(:,1:end-1) - repmat(h,Nvar,1) .* Y;
endfunction
function x = __bvp4c_solve__ (t, x, h, odefun, bcfun, Nvar, Nint, s)
fun = @( x ) ( [__bvp4c_fun_u__(t,
reshape(x(1:Nvar*(Nint+1)),Nvar,(Nint+1)),
reshape(x([1:Nvar*Nint*s]+Nvar*(Nint+1)),Nvar,Nint,s),
h,
s,
Nint,
Nvar)(:) ;
__bvp4c_fun_K__(t,
reshape(x(1:Nvar*(Nint+1)),Nvar,(Nint+1)),
reshape(x([1:Nvar*Nint*s]+Nvar*(Nint+1)),Nvar,Nint,s),
odefun,
h,
s,
Nint,
Nvar)(:);
bcfun(reshape(x(1:Nvar*(Nint+1)),Nvar,Nint+1)(:,1),
reshape(x(1:Nvar*(Nint+1)),Nvar,Nint+1)(:,end));
] );
x = fsolve ( fun, x );
endfunction
%!demo
%! a = 0;
%! b = 4;
%! Nint = 3;
%! Nvar = 2;
%! s = 3;
%! t = linspace(a,b,Nint+1);
%! h = diff(t);
%! u_1 = ones(1, Nint+1);
%! u_2 = 0*u_1;
%! u_0 = [u_1 ; u_2];
%! f = @(t,u) [ u(2); -abs(u(1)) ];
%! g = @(ya,yb) [ya(1); yb(1)+2];
%! solinit.x = t; solinit.y=u_0;
%! sol = bvp4cold(f,g,solinit);
%! [t,y]=ode45(f,[a b],[sol.y(:,1)]);
%! plot (sol.x,sol.y,'x-',t,y,'o-');
%!demo
%! a = 0;
%! b = 4;
%! Nint = 2;
%! Nvar = 2;
%! s = 3;
%! t = linspace(a,b,Nint+1);
%! h = diff(t);
%! u_1 = -ones(1, Nint+1);
%! u_2 = 0*u_1;
%! u_0 = [u_1 ; u_2];
%! f = @(t,u) [ u(2); -abs(u(1)) ];
%! g = @(ya,yb) [ya(1); yb(1)+2];
%! solinit.x = t; solinit.y=u_0;
%! sol = bvp4cold(f,g,solinit);
%! plot (sol.x,sol.y,'x-')