1 ## Copyright (C) 2008-2012 Carlo de Falco
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.
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.
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/>.
17 ## @deftypefn {Function File} {@var{A}} = bvp4c (@var{odefun}, @var{bcfun}, @var{solinit})
19 ## Solves the first order system of non-linear differential equations defined by
20 ## @var{odefun} with the boundary conditions defined in @var{bcfun}.
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:
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
30 ## @item @var{sol.solver} = "bvp4c" for compatibility
35 ## Author: Carlo de Falco <carlo@guglielmo.local>
36 ## Created: 2008-09-05
39 function sol = bvp4c(odefun,bcfun,solinit,options)
41 if (isfield(solinit,"x"))
44 error("bvp4c: missing initial mesh solinit.x");
47 if (isfield(solinit,"y"))
50 error("bvp4c: missing initial guess");
53 if (isfield(solinit,"parameters"))
54 error("bvp4c: solving for unknown parameters is not yet supported");
60 if (isfield(options,"RelTol"))
61 RelTol = options.RelTol;
63 if (isfield(options,"RelTol"))
64 AbsTol = options.AbsTol;
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);
84 du(:,kk) = odefun(t(kk), u(:,kk));
87 tm = (t(1:end-1)+t(2:end))/2;
90 um(nn,:) = interp1(t,u(nn,:),tm);
95 f_est(:,kk) = odefun(tm(kk), um(:,kk));
100 du_est(nn,:) = diff(u(nn,:))./h;
103 err = max(abs(f_est-du_est)); semilogy(tm,err), pause(.1)
105 RelErr = AbsErr/norm(du,inf)
107 if ( (AbsErr >= AbsTol) && (RelErr >= RelTol) )
108 ref_int = find( (err >= AbsTol) & (err./max(max(abs(du))) >= RelTol) );
113 t = sort([t, t_add]);
118 u_0(nn,:) = interp1(t_old, u(nn,:), t);
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);
139 sol.solver = 'bvp4c';
143 function diff_K = __bvp4c_fun_K__ (t, u, Kin, f, h, s, Nint, Nvar)
146 persistent C = [0 1/2 1 ];
148 persistent A = [0 0 0;
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);
163 function diff_u = __bvp4c_fun_u__ (t, u, K, h, s, Nint, Nvar)
166 persistent B= [1/6 2/3 1/6 ];
168 Y = zeros(Nvar, Nint);
170 Y += B(jj) * K(:,:,jj);
172 diff_u = u(:,2:end) - u(:,1:end-1) - repmat(h,Nvar,1) .* Y;
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),
185 reshape(x(1:Nvar*(Nint+1)),Nvar,(Nint+1)),
186 reshape(x([1:Nvar*Nint*s]+Nvar*(Nint+1)),Nvar,Nint,s),
192 bcfun(reshape(x(1:Nvar*(Nint+1)),Nvar,Nint+1)(:,1),
193 reshape(x(1:Nvar*(Nint+1)),Nvar,Nint+1)(:,end));
196 x = fsolve ( fun, x );
207 %! t = linspace(a,b,Nint+1);
209 %! u_1 = ones(1, Nint+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-')
224 %! t = linspace(a,b,Nint+1);
226 %! u_1 = -ones(1, Nint+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-')