ex9.m

fokker-planck example

timestep(1)=.1;
timestep(2)=.2;
timestep(3)=.4;
timestep(4)=.8;

mean=zeros(4,1);
eweaknorm=zeros(4,1);
meanew1000=0;

nsample=1000;

for j=1:nsample
 
  dtex=.1;
  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
  meanew1000=meanew1000+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  
  
  mean(1)=mean(1)+x(1000);


  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
 
  mean(2)=mean(2)+x(500);
  

  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  
  
  mean(3)=mean(3)+x(250);


  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  
  
  mean(4)=mean(4)+x(125);

end

eweaknorm=abs(mean/nsample-meanew1000/nsample); 

plot(timestep,eweaknorm,'-');
title('Weak Error versus Timestep  for Euler Scheme')
xlabel('Timestep Delta t')
ylabel('Absolute Weak (mean) Error at T=10')