Polyfit error warning for SDE

6 visualizzazioni (ultimi 30 giorni)
Mel Hernandez
Mel Hernandez il 5 Apr 2019
Commentato: Mel Hernandez il 5 Apr 2019
I'm getting this Warning everytime I use polyfit, and I can't seem to find the error behind my code. Here's the warning:
Warning: Polynomial is not unique; degree >= number of data points.
I'm trying to find the slope of 4 points from using the Euler-Maruyama method applied to my SDE. Here's the code:
% Euler-Maruyama method on linear SDE
% SDE is dX = -cos^3(X)*sin(X)*dt + cos^2(X)*dW
% X(0) = Xzero, Xzero = 0
% For 0 <= t <= 1
% Time steps: 1/50, 1/100, 1/150, 1/200
% Exact Solution: X(t) = tan^-1(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 = 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)*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 EM update & these increments
for j = 1:N_s
for k = 1:K
Xzero = 0;
X_app(j,1) = Xzero;
X_app(j,k+1) = X_app(j,k) + (-cos(X_app(j,k))^3)*(sin(X_app(j,k)))*tau + (cos(X_app(j,k))^2)*dW(j,k);
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) = atan(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 = 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)
I don't understand why I'm getting the error. It should be a 1 term polynomial (y = mx+b), and I only need the slope for my strong order of convergence. I get an answer with the warning, but I don't want to feel contempt with my answer when the warning is happening.

Risposta accettata

Jan
Jan il 5 Apr 2019
Modificato: Jan il 5 Apr 2019
tau and GErr are scalars, so you have one single point for the polynomial only. Then you cannot determine a slope.
I guess you want to create vectors instead. Then replace all tau and GErr to tau(f) and GErr(f):
...
tau(f) = T/K;
...
W(q,p+1) = W(q,p) + sqrt(tau(f))*normrnd(0,1);
...
GErr(f) = max(Err);
...
Then this will be working:
di = polyfit(log(tau),log(GErr),1);
By the way, this is not useful:
for k = 1:K
Xzero = 0; % ??? Set this once outside the loops
X_app(j,1) = Xzero; % ??? Why setting this in each iteration again
X_app(j,k+1) = X_app(j,k) + (-cos(X_app(j,k))^3)*(sin(X_app(j,k)))*tau + (cos(X_app(j,k))^2)*dW(j,k);
end
This:
for c = 1:N_s
for n = 1:K
dW(c,n) = W(c,n+1) - W(c,n);
end
end
can be simpified to:
dW = diff(W, 1, 2)
and
for tk = 1:K+1
for g = 1:N_s
X_exact(g,tk) = atan(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
to
X_exact = atan(W);
Xerr = abs(X_exact - X_app);
Err = sum(Xerr, 1) / N_s;
  1 Commento
Mel Hernandez
Mel Hernandez il 5 Apr 2019
Thank you so much! I completely forgot about tau(f) and GErr(f). It fixed it! And you didn't have to clean up my code. I do appreciate it though. Realized I did have a few too many lines going on. Thank you again!

Accedi per commentare.

Più risposte (0)

Categorie

Scopri di più su Time Series Objects in Help Center e File Exchange

Prodotti

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by