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?

score 2 · Answer 1 · answered Mar 24 '21 at 08:19

2

Your ExactY is a 1x5 array, whereas t is a 1x6 array (as is Y). For plot the input arrays for each x/y pair must have the same length, as indicated by the error message.

This line creates a 1x6 ExactY but it is redundant since ExactY is overwritten on the next line

ExactY = zeros(1, b+1);

It's unclear to me what the correct indices are, but something like this would resolve the error, perhaps only you can decide what the correct approach is but the same principal will apply:

ExactY=zeros(1,b+1);
ExactY(2:end)=(2^(1/3)+(1/3)*W).^3;

answered Mar 24 '21 at 08:19

Wolfie

27,562
7
28
55

@Julian as I said in my answer, I have no context from your question for what `ExactY` is *meant* to be! The error is your uneven array lengths, the reason for your uneven array lengths is because `W` and `t` are not the same length, and `W` generates `ExactY` - you need to review your maths to find where the issue is – Wolfie Mar 24 '21 at 08:32
Yeah you can use `t(1:5)` or `t(2:end)` to plot the 2nd line but that's very similar to how I padded `ExactY` - I have no clue whether it represents the correct answer. Note that you can't make `t` a `1x5` array completely because you're also using it to plot against the `1x6` array `Y` – Wolfie Mar 24 '21 at 09:28

Lutz Lehmann · Accepted Answer · 2021-03-26T14:56:54.353

2

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')

edited Mar 26 '21 at 14:56

answered Mar 24 '21 at 11:12

Lutz Lehmann

25,219
2
22
51

1

It is not well-seen if a question is modified in a way that the topic that is discussed disappears. Reactions to comments are, I believe, free game, but answers should remain anwers to the question. Better add a new section for the improved code. /// That said, you now have to remove `sd` and `deltaw` from the `eulermaruyama_solver`, and yes, if you shift the square root to the correct construction of the Brownian motion, then all other occurrences have to be removed. – Lutz Lehmann Mar 24 '21 at 11:41
1

@Julian you are treating this like a forum when it is not one, please resist editing questions to introduce scope creep - you asked a question about a specific error which is good, you've since transitioned this question to be about an issue with your algorithm which isn't ideal, now you're editing the original code which can make existing answers obsolete. Please read [ask], make edits clear to new visitors, and keep your questions specific – Wolfie Mar 24 '21 at 11:41
Yes, the only difference is a shift the place of operations. – Lutz Lehmann Mar 31 '21 at 14:39

Inconsistent vector length when using the `plot` command

2 Answers2