ex10.m

fokker-planck example

clear all;
timestep(1)=.01;
timestep(2)=.02;
timestep(3)=.04;
timestep(4)=.08;

estrongnorm=zeros(4,1);
mstrongnorm=zeros(4,1);

nsample=1;

for j=1:nsample
 
  dtex=.01;
  b=0.2;

  r=randn(1000,1);

  w0=0;
  w(1)=w0;

  for i=2:1000
    w(i)=w(i-1)+sqrt(dtex)*r(i);
  end  


  %the exact solution
  ew1000=exp(b*w(1000)); 

 
  %the numerical solution with the Euler scheme
  dt=timestep(1);
  x0=1.0;
  deltaw=sqrt(dtex)*r(1);
  x(1)=x0;
  for i=2:1000
    deltaw=sqrt(dtex)*r(i);
    x(i)=x(i-1)+0.5*b*b*x(i-1)*dt+x(i-1)*b*deltaw;
  end  
  
  estrongnorm(1)=estrongnorm(1)+abs(x(1000)-ew1000);
  
 %the numerical solution with the Milstein scheme
  dt=timestep(1);
  x0=1.0;
  deltaw=sqrt(dtex)*r(1);
  x(1)=x0;
  for i=2:1000
    deltaw=sqrt(dtex)*r(i);
    x(i)=x(i-1)+0.5*b*b*x(i-1)*dt+x(i-1)*b*deltaw+...
    0.5*b*b*x(i-1)*(deltaw*deltaw-dt);
  end 
  
  mstrongnorm(1)=mstrongnorm(1)+abs(x(1000)-ew1000);

  %the numerical solution with the Euler scheme
  dt=timestep(2);
  x0=1.0;
  deltaw=sqrt(dtex)*(r(1)+r(2));
  x(1)=x0;
  for i=2:500
    deltaw=sqrt(dtex)*(r(2*i-1)+r(2*i));
    x(i)=x(i-1)+0.5*b*b*x(i-1)*dt+x(i-1)*b*deltaw;
  end
  
  estrongnorm(2)=estrongnorm(2)+abs(x(500)-ew1000);
    
  %the numerical solution with the Milstein scheme
  dt=timestep(2);
  x0=1.0;
  deltaw=sqrt(dtex)*(r(1)+r(2));
  x(1)=x0;
  for i=2:500
    deltaw=sqrt(dtex)*(r(2*i-1)+r(2*i));
    x(i)=x(i-1)+0.5*b*b*x(i-1)*dt+x(i-1)*b*deltaw+...
    0.5*b*b*x(i-1)*(deltaw*deltaw-dt);
  end  

  mstrongnorm(2)=mstrongnorm(2)+abs(x(500)-ew1000);

  %the numerical solution with the Euler scheme
  dt=timestep(3);
  x0=1.0;
  deltaw=sqrt(dtex)*(r(1)+r(2)+r(3)+r(4));
  x(1)=x0;
  for i=2:250
    deltaw=sqrt(dtex)*(r(4*i-3)+r(4*i-2)+r(4*i-1)+r(4*i));
    x(i)=x(i-1)+0.5*b*b*x(i-1)*dt...
         +x(i-1)*b*deltaw;
  end  
  
  estrongnorm(3)=estrongnorm(3)+abs(x(250)-ew1000);

 %the numerical solution with the Milstein scheme
  dt=timestep(3);
  x0=1.0;
  deltaw=sqrt(dtex)*(r(1)+r(2)+r(3)+r(4));
  x(1)=x0;
  for i=2:250
    deltaw=sqrt(dtex)*(r(4*i-3)+r(4*i-2)+r(4*i-1)+r(4*i));
    x(i)=x(i-1)+0.5*b*b*x(i-1)*dt...
         +x(i-1)*b*deltaw+...
	 0.5*b*b*x(i-1)*(deltaw*deltaw-dt);
  end  
   
  mstrongnorm(3)=mstrongnorm(3)+abs(x(250)-ew1000);

  %the numerical solution with the Euler scheme
  dt=timestep(4);
  x0=1.0;
  deltaw=sqrt(dtex)*(r(1)+r(2)+r(3)+r(4)+r(5)+r(6)+r(7)+r(8));
  x(1)=x0;
  for i=2:125
    deltaw=sqrt(dtex)*(r(8*i-7)+r(8*i-6)+r(8*i-5)+r(8*i-4)+r(8*i-3)...
	 +r(8*i-2)+r(8*i-1)+r(8*i));
    x(i)=x(i-1)+0.5*b*b*x(i-1)*dt...
         +x(i-1)*b*deltaw;
  end  
  
  estrongnorm(4)=estrongnorm(4)+abs(x(125)-ew1000);
  
 %the numerical solution with the Milstein scheme
  dt=timestep(4);
  x0=1.0;
  deltaw=sqrt(dtex)*(r(1)+r(2)+r(3)+r(4)+r(5)+r(6)+r(7)+r(8));
  x(1)=x0;
  for i=2:125
    deltaw=sqrt(dtex)*(r(8*i-7)+r(8*i-6)+r(8*i-5)+r(8*i-4)+r(8*i-3)...
	 +r(8*i-2)+r(8*i-1)+r(8*i));
    x(i)=x(i-1)+0.5*b*b*x(i-1)*dt...
         +x(i-1)*b*deltaw+...
	 0.5*b*b*x(i-1)*(deltaw*deltaw-dt);
  end 
  
  mstrongnorm(4)=mstrongnorm(4)+abs(x(125)-ew1000);

end
 
plot(timestep,estrongnorm/nsample,'.',timestep,mstrongnorm/nsample,'-');
title('Strong Error versus Timestep for Euler & Milstein Schemes')
xlabel('Timestep Delta t')
ylabel('Absolute Stromng Error at T=10')
legend('Euler Scheme','Milstein Scheme')