I am trying to apply the Ito Taylor method to the SDE. I have attached what I'm using for the Ito Taylor because I would like for it to be double checked please.
I'm supposed to get a strong speed of convergence of 1.5, but I am getting a negative number, and I don't understand why.
SDE: dX = Xdt + XdW, X(0)=1
Exact: X(t) = exp(0.5*t + W(t))
I believe the problem lies in either my X_app and/or X_exact. Here's my code:
% Ito Taylor method on linear SDE
% SDE is dX = X*dt + X*dW
% X(0) = 1
% For 0 <= t <= 1
% Time steps: 1/50, 1/100, 1/150, 1/200
% Exact Solution: X(t) = exp(0.5*t + W(t))
% K = 50, 100, 150, 200; % number of steps
% T = 1; % end of interval
% tau = T/K; % time step
% A for loop for the 4 timesteps (1/50, 1/100, 1/150, 1/200)
for f = 1:4
K = 50*f; % So we can get 50, 100, 150, 200
T = 1;
tau(f) = T/K;
N_s = 2000; % number of trajectories
W = zeros(N_s,K+1); % setting up discrete Brownian path
dW = zeros(N_s,K); % setting up Brownian increments
X_app = zeros(N_s,K+1); % setting up Approx solution
X_exact = zeros(N_s,K+1); % setting up Exact solution
Xerr = zeros(1,K+1); % setting up Global error
% Compute 2000 independent Brownian Motion trajectories
for q = 1:N_s
for p = 1:K
W(q,p+1) = W(q,p) + sqrt(tau(f))*normrnd(0,1);
end
end
% Use these Brownian Motion trajectories to define the increments
for c = 1:N_s
for n = 1:K
dW(c,n) = W(c,n+1) - W(c,n);
end
end
% Generate 2000 solution trajectories using the IT update & these increments
for j = 1:N_s
for k = 1:K
X_app(j,1) = 1;
X_app(j,k+1) = X_app(j,k) + X_app(j,k)*tau(f) + X_app(j,k)*dW(j,k) +
X_app(j,k)*((tau(f))^2)*0.5 +
X_app(j,k)*tau(f)*dW(j,k) + X_app(j,k)*(0.5*
((dW(j,k)^2)-tau(f))) + X_app(j,k)*
(((1/6)*dW(j,k)^3)-(tau(f)*dW(j,k)*0.5));
end
end
% Estimate the global error of tau for the timestep using the exact
% solution
for tk = 1:K+1
for g = 1:N_s
X_exact(g,tk) = exp(0.5.*tk + W(g,tk));
Xerr(g,tk) = abs(X_exact(g,tk) - X_app(g,tk));
end
Err(tk) = (1/N_s)*sum(Xerr(:,tk));
end
GErr(f) = max(Err)
% Make a log-log plot of the data through a scatter plot
scatter(log(tau),log(GErr),'filled'), hold on
end
hold off
% Finding the slop of the line is the strong order of convergence
di = polyfit(log(tau),log(GErr),1);
slope = di(1)