Estimating Forcing Functions for the FitzHugh-Nagumo Equations
This page provides a demonstration of the use of forcing functiona perturbed set of FitzHugh-Nagumo Equations:
Here we will assume knowledge of $a$, $b$ and $c$, and estimate $g(t)$ from data.
The format of this demonstration follows that detailed in FhNEx.html and commentary will therefore be restricted to terms unique to forcing function estimation.
Contents
- RHS Functions
- Observation times
- Various Parameters
- Extra Information for the System:
- Penalties on Forcing Functions
- Initial Forcing Estimates
- Create trajectories
- Set up observations
- Fitting parameters
- Optimisation control
- Setting up functional data objects
- Smooth the data
- Re-Smooth with a Model-Based Penalty
- Now do the profiled estimation
- Examine Forcing Functions
- Calculate a Sample Information Matrix
- Calculate Hotelling Distance from Zero
RHS Functions
Since we are using a linear function to begin with, we make use of the forcing set of functions (although fhnfunodep will be used to produce data).
odefn = @fhnfunodep; % Function for ODE solver (exact) fn.fn = @forcingfun; % RHS function fn.dfdx = @forcingdfdx; % Derivative wrt inputs (Jacobian) fn.dfdp = @forcingdfdp; % Derviative wrt parameters fn.d2fdx2 = @forcingd2fdx2; % Hessian wrt inputs fn.d2fdxdp = @forcingd2fdxdp; % Cross derivatives wrt inputs and parameters fn.d2fdp2 = @forcingd2fdp2; % Hessian wrt parameters fn.d3fdx2dp = @forcingd3fdx2dp; % Third derivative wrt inputs, inputs, pars fn.d3fdx3 = @forcingd3fdx3; % Third derivative wrt inputs fn.d3fdxdp2 = @forcingd3fdxdp2; % Third derivative wrt inputs, pars and pars
Observation times
tspan = 0:0.05:20; % Observation times obs_pts{1} = 1:length(tspan); % Which components are observed at obs_pts{2} = 1:length(tspan); % which observation times. tfine = 0:0.05:20; % Times to plot solutions
Various Parameters
y0 = [-1,1]; % Initial conditions pars0 = [0.2; 0.2; 3]; % Parameters for the FitzHugh-Nagumo equations % Set up a perturbation functional data object basis_obj = create_bspline_basis([0 20],42,3,0:0.5:20); quadvals = MakeQuadPoints(0:0.5:20,21); % We will need to use basis_obj = putquadvals(basis_obj,quadvals); % quadrature points later pcoef = zeros(42,1); pcoef(floor(42/3):(ceil(42/3)+5)) = 10; fd_obj = fd(pcoef,basis_obj); % Finally decide on a noise level: sigma = 0.25; % Noise Level
Extra Information for the System:
In particular, we need to specify the original ODEs and their derifatives with respect to components (in this case the FitzHugh-Nagumo Equations) and also the option of adding further information for the original system.
fn_extras.fn = @fhnfun; % Original function fn_extras.dfdx = @fhndfdx; % First derivative fn_extras.d2fdx2 = @fhnd2fdx2; % Second derivative fn_extras.d3fdx3 = @fhnd3fdx3; % Third derivative fn_extras.extras = []; % Original information to fn_extras.fn. % We also need to add parameters: fn_extras.pars = pars0; % We also need to specify a basis for the forcing components and a vector % indicating which components are to be forced. fn_extras.basisp = {basis_obj}; % Cell array of basis functions. fn_extras.which = 1; % Force the first component of the % system only. % Note that although we are only forcing one component here, the forcing % basis is still represented as a cell array and we could equally have % estimated a number of forcing functions.
Penalties on Forcing Functions
We can also place roughness penalties on forcing functions, these will then occur as inputs into Profile_GausNewt.
pen = @forcingpen; % Penalty dpen = @forcingdpen; % Derivative with respect to forcing co-efficients d2pen = @forcingd2pen; % Second derivative % These penalty functions also require extra arguments in the form of a % struct to specify the basis, the degree of smoothing and the smoothing % penalty: pen_extras. basis = {basis_obj}; % Same basis as fn_extras. pen_extras.deg = 2; % Penalize the second derivative pen_extras.lambda = 0.0001; % Smoothing parameter.
Initial Forcing Estimates
We start off assuming the forcing function is zero. Since the coefficients of a basis expansion for it occur as parameters in the profiled estimation scheme we set them as:
startpars = zeros(getnbasis(basis_obj),1);
Create trajectories
odeopts = odeset('RelTol',1e-13);
[full_time,full_path] = ode45(odefn,tspan,y0,odeopts,pars0,fd_obj);
[plot_time,plot_path] = ode45(odefn,tfine,y0,odeopts,pars0,fd_obj);
Set up observations
Tcell = cell(1,size(full_path,2)); path = Tcell; for(i = 1:length(obs_pts)) Tcell{i} = full_time(obs_pts{i}); path{i} = full_path(obs_pts{i},i); end % add noise Ycell = path; for(i = 1:length(path)) Ycell{i} = path{i} + sigma*randn(size(path{i})); end % and set wts wts = []; if(isempty(wts)) % estimate wts if not given for(i = 1:length(Ycell)) wts(i) = 1./sqrt(var(Ycell{i})); end end
Fitting parameters
lambda = 1000; % Smoothing for model-based penalty lambda = lambda*wts; lambda0 = 1; % Smoothing for 1st-derivative penalty nknots = 401; % Number of knots to use. nquad = 5; % No. between-knots quadrature points. norder = 6; % Order of B-spline approximation
Optimisation control
lsopts_out = optimset('DerivativeCheck','off','Jacobian','on',... 'Display','iter','MaxIter',1000,'TolFun',1e-8,'TolX',1e-10); % Other observed optimiation control lsopts_other = optimset('DerivativeCheck','off','Jacobian','on',... 'Display','on','MaxIter',1000,'TolFun',1e-14,'TolX',1e-14,... 'JacobMult',@SparseJMfun); % Optimiation control within profiling lsopts_in = optimset('DerivativeCheck','off','Jacobian','on',... 'Display','off','MaxIter',1000,'TolFun',1e-14,'TolX',1e-14,... 'JacobMult',@SparseJMfun);
Setting up functional data objects
% set up knots range = [min(full_time),max(full_time)]; % range of observations knots_cell = cell(size(path)); % knots for each basis knots_cell(:) = {linspace(range(1),range(2),nknots)}; % set up bases basis_cell = cell(1,length(path)); % Create cell arrays. Lfd_cell = cell(1,length(path)); nbasis = zeros(length(path),1); bigknots = knots_cell{1}; % bigknots used for quadrature points nbasis(1) = length(knots_cell{1}) + norder - 2; for(i = 2:length(path)) bigknots = [bigknots knots_cell{i}]; nbasis(i) = length(knots_cell{i}) + norder -2; end quadvals = MakeQuadPoints(bigknots,nquad); % Create simpson's rule % quadrature points and values for(i = 1:length(path)) basis_cell{i} = MakeBasis(range,nbasis(i),norder,... % create bases knots_cell{i},quadvals,4); % with quadrature Lfd_cell{i} = fdPar(basis_cell{i},1,lambda0); % pts attatched end
Smooth the data
DEfd = smoothfd_cell(Ycell,Tcell,Lfd_cell); coefs = getcellcoefs(DEfd); devals = eval_fdcell(tfine,DEfd,0); for(i = 1:length(path)) subplot(length(path),1,i); plot(plot_time,plot_path(:,i),'b','LineWidth',2); hold on; plot(tfine,devals{i},'r','LineWidth',2); plot(Tcell{i},Ycell{i},'b.'); hold off; end
Re-Smooth with a Model-Based Penalty
[newcoefs,resnorm2] = lsqnonlin(@SplineCoefErr,coefs,[],[],... lsopts_other,basis_cell,Ycell,Tcell,wts,lambda,fn,[],startpars,... fn_extras); tDEfd = Make_fdcell(newcoefs,basis_cell); % plot results along with exact solution, there is a noticeable lack of % fit. devals = eval_fdcell(tfine,tDEfd,0); for(i = 1:length(Ycell)) subplot(length(Ycell),1,i); plot(tfine,devals{i},'r','LineWidth',2); hold on; plot(Tcell{i},Ycell{i},'b.'); plot(plot_time,plot_path(:,i),'c'); hold off end
Optimization terminated: relative function value changing by less than OPTIONS.TolFun.
Now do the profiled estimation
Recall that at this point we are estimating the coefficients of the forcing functions.
[fcoefs,DEfd] = Profile_GausNewt(startpars,lsopts_out,tDEfd,fn,lambda,... Ycell,Tcell,wts,[],lsopts_in,fn_extras,pen,dpen,pen_extras); % DEfd now takes the form of a smooth to the data and we can see that it % fits better: devals = eval_fdcell(tfine,DEfd,0); for(i = 1:length(path)) subplot(length(path),1,i); plot(tfine,devals{i},'r','LineWidth',2); hold on; plot(Tcell{i},Ycell{i},'b.'); plot(plot_time,plot_path(:,i),'c'); hold off end
Optimization terminated: relative function value
changing by less than OPTIONS.TolFun.
Iteration steps Residual Improvement Grad-norm parameters
1 1 125.534 0.572156 0.375 -2.1382 -0.12759 1.0913 -0.76994 -0.94532 -0.3426 0.15708 0.65152 -0.98157 -0.9668 -0.92148 -0.55023 -0.48966 7.8561 -2.9909 1.4163 2.5393 4.269 8.9936 -2.4684 9.1373 7.1025 -1.4295 0.93904 0.22695 1.5873 0.49971 -0.19859 -1.0754 2.2416 -2.4401 -4.3036 2.1472 -0.66319 0.45863 -0.074892 -0.11846 0.25127 -0.064283 0.24839 -0.40145 0.19478
2 1 72.9574 0.418825 0.0641 -1.79122 0.130561 1.14831 -0.992914 -1.11579 -0.59249 -0.210706 0.311705 -1.28529 -0.284699 -0.462374 0.850836 -0.350049 15.0343 1.61657 0.546123 6.29447 10.5022 9.17969 -2.19989 7.57749 7.95961 0.442467 -0.539036 -0.20414 1.08514 0.468304 0.0964572 0.0526961 0.033122 0.11396 -2.983 1.32087 -0.0861879 0.525705 0.104564 -0.0312305 0.206797 -0.229667 0.249042 -0.155339 -0.00186128
3 1 45.0095 0.383071 0.0105 -2.17465 -0.131696 1.10738 -0.733679 -0.844844 -0.28101 0.166931 0.804544 -0.882374 -0.107278 -0.418309 0.134699 0.78076 7.06232 9.88573 3.21205 7.95139 12.0219 9.97113 -1.52484 5.20258 6.22403 -0.315997 -1.16938 -0.805447 0.784912 0.304572 -0.0178085 -0.511926 0.381613 0.593375 -2.40877 1.54383 -0.337496 0.504443 0.0620958 -0.0346536 0.231071 -0.143 0.288254 -0.224664 -0.0372753
4 1 39.033 0.132783 0.00103 -2.1732 -0.129636 1.10649 -0.721991 -0.836516 -0.260538 0.193023 0.826744 -0.765653 -0.0344533 -0.357258 0.291074 0.329352 8.81437 10.5126 8.9299 9.89112 10.9539 10.3127 -0.859392 2.31127 2.65512 -0.697649 -0.682689 -0.602323 0.948964 0.419219 0.0577328 -0.581469 0.258023 0.467999 -1.79752 1.5636 -0.310146 0.53662 0.0729478 -0.0384608 0.232991 -0.148714 0.279947 -0.249074 -0.0548577
5 1 38.732 0.00771044 0.0002 -2.17756 -0.131109 1.10605 -0.724934 -0.837979 -0.263995 0.189652 0.819924 -0.76845 -0.0426131 -0.356336 0.281817 0.35209 8.76981 10.4791 10.3263 10.6891 10.748 10.2451 -0.564259 1.09551 1.26145 -0.83774 -0.598384 -0.638294 0.946939 0.412166 0.0494583 -0.587241 0.251499 0.478445 -1.77983 1.55306 -0.305222 0.534926 0.0725833 -0.0370508 0.230586 -0.146889 0.277929 -0.245776 -0.0549932
6 1 38.7302 4.75136e-005 1.05e-006 -2.1776 -0.131137 1.10602 -0.72512 -0.838121 -0.264236 0.189344 0.819722 -0.769168 -0.0424373 -0.357133 0.282624 0.351 8.77238 10.4872 10.3622 10.7155 10.7493 10.2418 -0.519915 0.914609 1.18221 -0.839014 -0.592998 -0.639154 0.946173 0.411611 0.047248 -0.584712 0.254381 0.468125 -1.76516 1.54827 -0.305555 0.534735 0.0726157 -0.0377534 0.231464 -0.147232 0.279119 -0.247039 -0.0545797
7 1 38.7302 1.93535e-007 1.08e-007 -2.1776 -0.131136 1.10603 -0.725117 -0.838113 -0.264229 0.189345 0.819758 -0.769193 -0.0424019 -0.357144 0.282618 0.351088 8.77206 10.4877 10.3623 10.7156 10.7497 10.2412 -0.515876 0.899969 1.18327 -0.838322 -0.59293 -0.639203 0.946124 0.411487 0.0472336 -0.584726 0.254789 0.467344 -1.76436 1.5478 -0.305567 0.534761 0.072628 -0.0377672 0.231473 -0.147239 0.279153 -0.247097 -0.0545941
8 1 38.7302 9.34245e-010 1.27e-009 -2.1776 -0.131136 1.10603 -0.725116 -0.838113 -0.264229 0.189345 0.819759 -0.769193 -0.0424015 -0.357143 0.282618 0.351086 8.77206 10.4878 10.3623 10.7156 10.7497 10.2412 -0.515599 0.898941 1.18339 -0.838292 -0.592928 -0.6392 0.946124 0.41148 0.0472321 -0.584728 0.254811 0.467306 -1.76432 1.54779 -0.305573 0.534765 0.0726288 -0.0377691 0.231475 -0.14724 0.279156 -0.247103 -0.0545952
Examine Forcing Functions
First we constitute the approximated forcing function:
fd_approx = fd(fcoefs,basis_obj); % Now plot the results along with the original forcing function: force_devals = eval_fd(tfine,fd_obj); approx_devals = eval_fd(tfine,fd_approx); subplot(1,1,1) plot(tfine,force_devals,'r','LineWidth',2); hold on plot(tfine,approx_devals,'b','LineWidth',2); hold off
Calculate a Sample Information Matrix
% We know this will be wrong, but just ot see what this looks like. d2Jdp2 = make_d2jdp2(DEfd,fn,Ycell,Tcell,lambda,fcoefs,[],wts,... fn_extras,d2pen,pen_extras); d2JdpdY = make_d2jdpdy(DEfd,fn,Ycell,Tcell,lambda,fcoefs,[],wts,fn_extras); dpdY = -d2Jdp2\d2JdpdY; S = make_sigma(DEfd,Tcell,Ycell,0); Cov = dpdY*S*dpdY';
Calculate Hotelling Distance from Zero
Look at the distance of fcoefd from 0 with respect to the metric defined by Cov. This is a heuristic test for the goodness of fit of the original equations.
fcoefs'*inv(Cov)*fcoefs
ans = 1.8282e+003