function fig8_17

%  fit solution of IVP to data to find parameters p1 and p2
%  requires file splineA.m

%  IVP:   y'=ff(y,v)  v'=gg(y,v)  for 0 < t < tmax
%         y(0)=y0     v(0)=v0

%  method: uses RK4 + cubic spline + fminsearch
%  as written, the FitzHugh-Nagumo equations are used
%  the functions ff(y,v) and gg(y,v) are entered in function f(t,y,p1,p2)

global td zd tmax

%%%%  generate data for regression
%%%%  note that tmax, y0, v0 must also be entered in function F
%  first find numerical solution for given values
p1=3; p2=0.2; p3=0.2;
tmax=20;
y0=-1;
v0=1;
nt=100;
t=linspace(0,tmax,nt);
[t,z] = ode45(@(t,y) f(t,y,p1,p2,p3),t,[y0 v0]);
dz0=f(0,z(1,:),p1,p2,p3);
dz1=f(tmax,z(nt,:),p1,p2,p3);
% now randomize data:  nd = number of data points
% nd=120;
nd=60;
yd=zeros(1,nd); vd=zeros(1,nd); zd=zeros(nd,2);
%td=linspace(0.1,19.9,nd);
tdd=tmax*rand(1,nd);
td=sort(tdd);
td(1)=0;
qd=0.4*(rand(1,nd)-0.5);
qdd=0.5*(rand(1,nd)-0.5);
size(t)
size(z(:,1)')
size(td)
yd = splineA(t',z(:,1)',td,[1 1],[dz0(1) dz1(1)]);
vd = splineA(t',z(:,2)',td,[1 1],[dz0(2) dz1(2)]);
for in=2:nd
    yd(in)=(1+qd(in))*yd(in);
    vdd=vd(in);
    if abs(qdd(in))>0.2
        vd(in)=(1+qdd(in))*vd(in);
    else
        vd(in)=vd(in)+0.4*(rand-0.5);
    end
    zd(in,1)=yd(in);
    zd(in,2)=vd(in);
end

[sol,fv,ee,o] = fminsearch(@F, [1, 1, 1],optimset('MaxFunEvals',1000000,'MaxIter',1000000,'TolX',1e-6))
pd1=sol(1); pd2=sol(2); pd3=sol(3);

%  p1=c  p2=b  p3=a
% z(1)=p1*(y(1)-y(1)^3+y(2));
% z(2)=-(y(1)-p3+p2*y(2))/p1;
fprintf('\n  Computed Solution =  %e  %e   %e\n\n',pd1, pd2, pd3);

%%% plot data and solution using calcuated parameter values
clf
% get(gcf)
set(gcf,'Position', [4 1115 658 230])
subaxis(1,1,1,1,'MT',0.01,'MB',0.16,'MR',-0.02,'ML',0.05,'P',0.04)
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
nt=200;
t=linspace(0,tmax,nt);
[t,z] = ode45(@(t,y) f(t,y,pd1,pd2,pd3),t,[y0 v0]);
plot(t,z(:,1),'r','LineWidth',1.6);
hold on
plot(t,z(:,2),'b','LineWidth',1.6)
axis([0 tmax -1.5 1.5])
legend(' v',' w','Location','SouthWest','AutoUpdate','off','FontSize',16,'FontWeight','bold')
grid on
set(gca,'ytick',[-1.5 -0.5 0.5 1.5])
plot(td,yd,'.r','MarkerSize',18)
plot(td,vd,'db','MarkerSize',5,'LineWidth',1.5)
xlabel('t-axis')
%ylabel('y-axis')
set(gca,'FontSize',16,'FontWeight','bold')

%exportgraphics(gca,'/Users/mark/Desktop/fnff.eps')


%%% error function: F = (y-yd)^2+(v-vd)^2
% IMPORTANT:  m = number of points for cubic spline
function s=F(p)
global td zd tmax
p1=p(1); p2=p(2); p3=p(3);
m=60;
y0=-1;
v0=1;
z0=[y0 v0]';
t=linspace(0,tmax,m);
k=t(2)-t(1);
[y,v]=rk4(t,z0,k,m,p1,p2,p3);
dy0=f(0,z0,p1,p2,p3);
dy1=f(tmax,z0,p1,p2,p3);
yf = fspline(t,y,dy0,td);
vf = fspline(t,v,dy1,td);
s=norm(yf-zd(:,1)',2)^2+norm(vf-zd(:,2)',2)^2;

% right-hand side of DE
function z=f(t,y,p1,p2,p3)
z=zeros(2,1);
z(1)=p1*(y(1)-y(1)^3+y(2));
z(2)=-(y(1)-p3+p2*y(2))/p1;

% RK4 method
function [ypoints, vpoints] =rk4(t,y0,h,n,p1,p2,p3)
y=y0;
ypoints=y0(1);
vpoints=y0(2);
for i=2:n
    k1=h*f(t(i-1),y,p1,p2,p3);
    k2=h*f(t(i-1)+0.5*h,y+0.5*k1,p1,p2,p3);
    k3=h*f(t(i-1)+0.5*h,y+0.5*k2,p1,p2,p3);
    k4=h*f(t(i),y+k3,p1,p2,p3);
    yy=y+(k1+2*k2+2*k3+k4)/6;
    ypoints=[ypoints, yy(1)];
    vpoints=[vpoints, yy(2)];
    y=yy;
end































