Inconsistent vector length when using the `plot` command

Question

I have created a function below. The function is a generic solver for a stochastic differential equation that uses the Euler-Maruyama method.

function [Y,t]=eulermaruyama_solver(Y0,T,a,b,c,F,G,deltaW)
% deltaW is called on in the test script

dt=T/b; t=[0:dt:T]; sd=sqrt(dt);
dw = zeros(a,b+1);
dw(:,2:b+1) = sd * deltaW;  
Y=zeros(c,b+1);
Y(:,1) = Y0; Yn = Y0;
for n = 1:b
    Y(:,n+1) = Yn + F(t(n),Yn)*dt + G(t(n),Yn) * dw(:,n+1);
    Yn = Y(:,n+1);
end

To test my function in the simple case that a=c=1, and I define functions F and G. Below is my test script.

a=1; c=1; Y0=1; T=1; b=5;
F=@(t,y) y^(1/2); G=@(t,y) 2*y;



% Euler-Maruyama solution
[deltaW,W]=randn_numbers(a,b);
[Y,t]=eulermaruyama_solver(Y0,T,a,b,c,F,G,deltaW);

% Exact solution
ExactY=zeros(1,b+1);
ExactY=(2^(1/5)+(1/2)*W).^3;

% Plot 
subplot(1,2,1)
plot(t, Y, 'r', t, ExactY,'b')

But the following error message is displayed: "Error using plot, Vectors must be the same length".

Have I made a mistake?

Answer 1

Mathematically, what goes wrong is that your exact solution is wrong on two counts.

• The most severe is that you do not apply the scaling with sqrt(dt), giving the wrong magnitudes.
• Less severe, but still giving a significant error is that your initial values do not match. Your exact solution Y(t)=(2^(1/3)+W(t)/3)^3 has Y(0)=2, however you use Y0=1. Perhaps you should set ExactY(2:b+1,:)=(Y0^(1/3)+(1/3)*W*sqrt(dt)).^3;, do not forget so set ExactY(1,:)=Y0;. Or construct W with one element more so that W(1,:)==0.

Btw., what do you suppose should happen in eulermaruyama_solver if a and c do not have the same value?

---

Keeping the structure, you could change the methods as

function [dW,W]=randn_numbers(a,b,T)
  W=zeros(a,b+1);
  dW=randn(a,b)*sqrt(T/b);
  W(:,2:end)=cumsum(dW,2);
end%function

function [Y,t]=eulermaruyama_solver(Y0,T,a,b,c,F,G,dW)
  dt=T/b; t=[0:dt:T]; 
  Y=zeros(a,b+1);
  Y(:,1) = Y0; Yn = Y0;
  for n = 1:b
    Y(:,n+1) = Yn + F(t(n),Yn)*dt + G(t(n),Yn) .* dW(:,n);
    Yn = Y(:,n+1);
  end%for
end%function

a=4; c=a; Y0=1; T=1; b=50;
F=@(t,y) (1/3)*y.^(1/3); G=@(t,y) y.^(2/3);

[dW,W]=randn_numbers(a,b,T);
[Y,t]=eulermaruyama_solver(Y0,T,a,b,c,F,G,dW);

% Exact solution
ExactY=(Y0^(1/3)+(1/3)*W).^3;

% Plot 
plot(t, Y, 'r', t, ExactY,'b')