octave_packages/odepkg-0.8.2/bvp4c.m

   1 ## Copyright (C) 2008-2012 Carlo de Falco
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 2 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {@var{A}} = bvp4c (@var{odefun}, @var{bcfun}, @var{solinit})
  18 ##
  19 ## Solves the first order system of non-linear differential equations defined by
  20 ## @var{odefun} with the boundary conditions defined in @var{bcfun}.
  21 ##
  22 ## The structure @var{solinit} defines the grid on which to compute the
  23 ## solution (@var{solinit.x}), and an initial guess for the solution (@var{solinit.y}).
  24 ## The output @var{sol} is also a structure with the following fields:
  25 ## @itemize
  26 ## @item  @var{sol.x} list of points where the solution is evaluated
  27 ## @item  @var{sol.y} solution evaluated at the points @var{sol.x}
  28 ## @item  @var{sol.yp} derivative of the solution evaluated at the
  29 ## points @var{sol.x}
  30 ## @item  @var{sol.solver} = "bvp4c" for compatibility
  31 ## @end itemize
  32 ## @seealso{odpkg}
  33 ## @end deftypefn
  34
  35 ## Author: Carlo de Falco <carlo@guglielmo.local>
  36 ## Created: 2008-09-05
  37
  38
  39 function sol = bvp4c(odefun,bcfun,solinit,options)
  40
  41   if (isfield(solinit,"x"))
  42     t = solinit.x;
  43   else
  44     error("bvp4c: missing initial mesh solinit.x");
  45   end
  46
  47   if (isfield(solinit,"y"))
  48     u_0 = solinit.y;
  49   else
  50     error("bvp4c: missing initial guess");
  51   end
  52
  53   if (isfield(solinit,"parameters"))
  54     error("bvp4c: solving for unknown parameters is not yet supported");
  55   end
  56
  57   RelTol = 1e-3;
  58   AbsTol = 1e-6;
  59   if ( nargin > 3 )
  60     if (isfield(options,"RelTol"))
  61       RelTol = options.RelTol;
  62     endif
  63     if (isfield(options,"RelTol"))
  64       AbsTol = options.AbsTol;
  65     endif
  66   endif
  67
  68   Nvar = rows(u_0);
  69   Nint = length(t)-1;
  70   s    = 3;
  71   h    = diff(t);
  72
  73   AbsErr  = inf;
  74   RelErr  = inf;
  75   MaxIt   = 10;
  76
  77   for iter = 1:MaxIt
  78
  79     x       = [ u_0(:); zeros(Nvar*Nint*s,1) ];
  80     x       = __bvp4c_solve__ (t, x, h, odefun, bcfun, Nvar, Nint, s);
  81     u       = reshape(x(1:Nvar*(Nint+1)),Nvar,Nint+1);
  82
  83     for kk=1:Nint+1
  84       du(:,kk) = odefun(t(kk), u(:,kk));
  85     end
  86
  87     tm = (t(1:end-1)+t(2:end))/2;
  88     um = [];
  89     for nn=1:Nvar
  90       um(nn,:) = interp1(t,u(nn,:),tm);
  91     endfor
  92
  93     f_est = [];
  94     for kk=1:Nint
  95       f_est(:,kk) = odefun(tm(kk), um(:,kk));
  96     end
  97
  98     du_est = [];
  99     for nn=1:Nvar
 100       du_est(nn,:) = diff(u(nn,:))./h;
 101     end
 102
 103     err    = max(abs(f_est-du_est)); semilogy(tm,err), pause(.1)
 104     AbsErr = max(err)
 105     RelErr = AbsErr/norm(du,inf)
 106
 107     if    ( (AbsErr >= AbsTol) && (RelErr >= RelTol) )
 108       ref_int = find( (err >= AbsTol) & (err./max(max(abs(du))) >= RelTol) );
 109
 110       t_add = tm(ref_int);
 111       t_old = t;
 112
 113       t     = sort([t, t_add]);
 114       h     = diff(t);
 115
 116       u_0 = [];
 117       for nn=1:Nvar
 118         u_0(nn,:) = interp1(t_old, u(nn,:), t);
 119       end
 120       Nvar = rows(u_0);
 121       Nint = length(t)-1
 122     else
 123       break
 124     end
 125
 126   endfor
 127
 128   ## K    = reshape(x([1:Nvar*Nint*s]+Nvar*(Nint+1)),Nvar,Nint,s);
 129   ## K1 = reshape(K(:,:,1), Nvar, Nint);
 130   ## K2 = reshape(K(:,:,2), Nvar, Nint);
 131   ## K3 = reshape(K(:,:,3), Nvar, Nint);
 132
 133
 134
 135   sol.x = t;
 136   sol.y = u;
 137   sol.yp= du;
 138   sol.parameters = [];
 139   sol.solver = 'bvp4c';
 140
 141 endfunction
 142
 143 function diff_K = __bvp4c_fun_K__ (t, u, Kin, f, h, s, Nint, Nvar)
 144
 145   %% coefficients
 146   persistent C = [0      1/2    1 ];
 147
 148   persistent A = [0      0      0;
 149                   5/24   1/3   -1/24;
 150                   1/6    2/3    1/6];
 151
 152   for jj = 1:s
 153     for kk = 1:Nint
 154       Y = repmat(u(:,kk),1,s) + ...
 155           (reshape(Kin(:,kk,:),Nvar,s) * A.') * h(kk);
 156       diff_K(:,kk,jj) = Kin(:,kk,jj) - f (t(kk)+C(jj)*h(kk), Y);
 157     endfor
 158   endfor
 159
 160 endfunction
 161
 162
 163 function diff_u = __bvp4c_fun_u__ (t, u, K, h, s, Nint, Nvar)
 164
 165   %% coefficients
 166   persistent B= [1/6 2/3 1/6 ];
 167
 168   Y = zeros(Nvar, Nint);
 169   for jj = 1:s
 170     Y +=  B(jj) * K(:,:,jj);
 171   endfor
 172   diff_u = u(:,2:end) - u(:,1:end-1) - repmat(h,Nvar,1) .* Y;
 173
 174 endfunction
 175
 176 function x = __bvp4c_solve__ (t, x, h, odefun, bcfun, Nvar, Nint, s)
 177   fun = @( x ) ( [__bvp4c_fun_u__(t,
 178                                   reshape(x(1:Nvar*(Nint+1)),Nvar,(Nint+1)),
 179                                   reshape(x([1:Nvar*Nint*s]+Nvar*(Nint+1)),Nvar,Nint,s),
 180                                   h,
 181                                   s,
 182                                   Nint,
 183                                   Nvar)(:) ;
 184                   __bvp4c_fun_K__(t,
 185                                   reshape(x(1:Nvar*(Nint+1)),Nvar,(Nint+1)),
 186                                   reshape(x([1:Nvar*Nint*s]+Nvar*(Nint+1)),Nvar,Nint,s),
 187                                   odefun,
 188                                   h,
 189                                   s,
 190                                   Nint,
 191                                   Nvar)(:);
 192                   bcfun(reshape(x(1:Nvar*(Nint+1)),Nvar,Nint+1)(:,1),
 193                         reshape(x(1:Nvar*(Nint+1)),Nvar,Nint+1)(:,end));
 194                   ] );
 195
 196   x    = fsolve ( fun, x );
 197 endfunction
 198
 199
 200
 201 %!demo
 202 %! a            = 0;
 203 %! b            = 4;
 204 %! Nint         = 3;
 205 %! Nvar         = 2;
 206 %! s            = 3;
 207 %! t            = linspace(a,b,Nint+1);
 208 %! h            = diff(t);
 209 %! u_1          = ones(1, Nint+1);
 210 %! u_2          = 0*u_1;
 211 %! u_0          = [u_1 ; u_2];
 212 %! f            = @(t,u) [ u(2); -abs(u(1)) ];
 213 %! g            = @(ya,yb) [ya(1); yb(1)+2];
 214 %! solinit.x = t; solinit.y=u_0;
 215 %! sol = bvp4c(f,g,solinit);
 216 %! plot (sol.x,sol.y,'x-')
 217
 218 %!demo
 219 %! a            = 0;
 220 %! b            = 4;
 221 %! Nint         = 2;
 222 %! Nvar         = 2;
 223 %! s            = 3;
 224 %! t            = linspace(a,b,Nint+1);
 225 %! h            = diff(t);
 226 %! u_1          = -ones(1, Nint+1);
 227 %! u_2          = 0*u_1;
 228 %! u_0          = [u_1 ; u_2];
 229 %! f            = @(t,u) [ u(2); -abs(u(1)) ];
 230 %! g            = @(ya,yb) [ya(1); yb(1)+2];
 231 %! solinit.x = t; solinit.y=u_0;
 232 %! sol = bvp4c(f,g,solinit);
 233 %! plot (sol.x,sol.y,'x-')