# What is the way of solving the noise added differential equations in matlab?

34 views (last 30 days)
Ms. Sita on 9 Feb 2022
Commented: Ms. Sita on 11 Feb 2022
I want to study the stoichastic resonance in duffing oscillator. The duffing oscillator is given as, ,
where are constant parameter and f is the forcing amplitude and is the white gaussian noise and D is the varience. I want to see the effect of noise on the dyanmics of the system. I tried solving the above equation using ode45 but it's not giving any results......
clear all;
clc;
clf;
%% ----------------INPUT PARAMETERS---------------------------
f=0.33; alpha=0.5; w0=sqrt(-1); beta=1; w=0.96; D=10;
x0 = [-1 0 0]; %initial condition
tspan = [0:0.1:100]; %time duration
n_rec= 300; %recording time
options = odeset ('RelTol',1e-7,'AbsTol',1e-6); %acuuracy set
[t,x] = ode45('duffing_sr',tspan,x0, options); %solving using ode45
[row, col]=size(x);
%Storing the data after discarding initial transients
x1 = x(n_rec:row,1);
x2 = x(n_rec:row,2);
x3 = x(n_rec:row,3);
n = length(x1); % Length of x1 after discarding the initial transients
t_n=t(n_rec:row);
plot(t_n,x1); %plotting
%---------------------------Function-----------------------------------------------
function dxdt = duffing_sr(t,x)
global f alpha w0 beta w D ;
dxdt = zeros(3,1);
dxdt(1) = x(2);
dxdt(2) = f*sin(x(3))-alpha*x(2)-w0^2*x(1)-beta*x(1)^3+sqrt(D)*(randn()-0.5);
dxdt(3) = w;
end
****************************************************************************************
I have no idea how to solve noise added differential equations.

Alan Stevens on 9 Feb 2022
Try to avoid the (mis)use of global!
%% ----------------INPUT PARAMETERS---------------------------
f=0.33; alpha=0.5; w0=sqrt(-1); beta=1; w=0.96; D=10;
x0 = [-1 0]; %initial condition
tspan = 0:0.1:100; %time duration
n_rec= 300; %recording time
%options = odeset ('RelTol',1e-7,'AbsTol',1e-6); %acuuracy set
[t,x] = ode45(@(t,x) duffing_sr(t,x,f,alpha,beta,w,D),tspan,x0); %solving using ode45
[row, col]=size(x);
%Storing the data after discarding initial transients
x1 = x(n_rec:row,1);
x2 = x(n_rec:row,2);
n = length(x1); % Length of x1 after discarding the initial transients
t_n=t(n_rec:row);
plot(t_n,x1); %plotting
%---------------------------Function---------------------------------------
function dxdt = duffing_sr(t,x,f,alpha,beta,w,D)
dxdt = [x(2);
f*sin(w*t)-alpha*x(2)-(-1)*x(1)-beta*x(1)^3+sqrt(D)*(randn-0.5)];
end
##### 2 CommentsShowHide 1 older comment
Ms. Sita on 11 Feb 2022
Thank you Alan Stevens for correcting the code.