A second order differential equation about brownian motion

I have two equations and (η is a random number, but η can't be repeatedly used. Otherwise, some movement tracks may be the same.). These are two equations about Brownian motion. Let's say there is only one particle. The first motion should be simulated for 1000 times, while the second one should be simulated for only once. How could I do? The result of both simulations needs to be shown in two histograms respectively. Also, I know that when solving this question, the equation of and x is needed. I want to know what exactly the two equations are.

Here's what I have done until now. Any help will be appreciated.

%%
m=1
r=1 % I am not sure whether these two suppositions are right or not
x(1)=1
v(1)=0
h=1
for i=(1:1000)
    a(i)=(r*v(i)+i)/m
    v(i)=v(1)+h*a(i)
    x(i)=x(1)+h*v(i)
end
hist(x(i),25)

%%
m=1
r=1 
k=1
h=1 % I am not sure whether these four suppositions are right or not
x(1)=1
v(1)=0
for i=(1:1000)
    a(i)=(r*v(i)-k*x(i)+i)/m
    v(i)=v(1)+h*a(i)
    x(i)=x(1)+h*v(i)
end
hist(x(i),25)

Answer 1

[T.Nikhil kumar]

Hello!

I understand from your question that you want to simulate the two equations of Brownian motion for one particle using MATLAB, them being:

ma = rv+η (to simulate for 1000 times)

ma = rv+kx+η (to simulate only once)

Your requirement is also to represent these equations in terms of dx/dt and x and ultimately plot these values of x obtained on separate histograms. I also understand that you want the noise term (η) to be a random number but not the same value for every time the simulation is run. Please find below a possible solution to your question.

Note:

1. We can rewrite the equations as:

m(dv/dt) = λ(dx/dt) + η

m(dv/dt) = λ(dx/dt) + kx + η

2. Here, we assume that the noise term (η) is independent of time.

3. We assume constant values for time step, dt=0.1 and Total simulation time (Stop time) as 10 s.

% Defining parameters and initial values of velocity 
r = 1;       % Damping coefficient (as mentioned in your code)
k = 1;       % Spring constant (as mentioned in your code)
m = 1;       % Mass of particle (as mentioned in your code)
v = 1;       % Initial velocity (as mentioned in your code)
dt = 0.1;    % Time step (step size for solving differential equation using Euler’s method /h )
T = 10;      % Total Simulation time

% Creating two variables denoting the number of times the equation must be simulated
numOfSimFirstEq=1000;
numOfSimSecEq=1;

% Creating two column vectors (initialized with 0’s as placeholders) 
% that will hold the results (positions of particle) of each simulation of the two equations
positionsFirstEq = zeros(numOfSimFirstEq, 1);
positionsSecondEq = zeros(numOfSimSecEq, 1);

% Simulating the first equation for 1000 times
for i = 1:numOfSimFirstEq

    % Generating normally distributed random values for noise term
    % (with a size of [1,1000] i.e unique for each simulation )
    noiseTerm = randn([1,1000]);
    
    x = 1;  % Initial position of the particle (as mentioned in your code)
    
    for t = 0:dt:T

        % For each time step, finding the position and velocity 
        % and updating it using the first equation (ma=rv+η)

        % find acceleration using a=(rv+η)/m
        a = (r * v + noiseTerm(i))/m;

        % update velocity by adding the change in velocity term to the
        % current velocity. It is known that a=dv/dt => dv=a.dt
        v = v + a * dt;
        
        % update position by adding the change in position term to the
        % current position. It is known that v=dx/dt => dx=v.dt
        x = x + v * dt;
    end
    
    % Storing the result of the first equation simulation
    positionsFirstEq(i) = x;
end

% Setting the value of Initial velocity again to 1
v=1;

% Simulating the second equation once
for i = 1:numOfSimSecEq

     % Generating normally distributed random values for noise term
    % (with a size of [1,1] i.e only for one simulation )
    noiseTerm = randn(1);
    
    x = 1;  % Initial position of the particle
    
    for t = 0:dt:T
       
        % For each time step finding the position and velocity 
        % and updating it using the first equation (ma=rv+η+kx)

        % find acceleration using a=(rv+η+kx)/m
        a = (r * v + noiseTerm + k * x)/m;

        % update velocity by adding the change in velocity term to the
        % current velocity. It is known that a=dv/dt => dv=a.dt
        v = v + a * dt;
        
        % update position by adding the change in position term to the
        % current position. It is known that v=dx/dt => dx=v.dt
        x = x + v * dt;
    end
    
    % Storing the result of the second equation simulation
    positionsSecondEq(i) = x;
end

% Plotting histograms side by side using subplot function
figure;

subplot(2, 1, 1); 

histogram(positionsFirstEq,25); % Number of bins set to 25 as mentioned in your code
title('Histogram for First Equation Simulation');
xlabel('Position');
ylabel('Probability Density');

subplot(2, 1, 2);

histogram(positionsSecondEq,25);
title('Histogram for Second Equation Simulation');
xlabel('Position');
ylabel('Probability Density');

Refer to the attached question about solving differential equations via Euler’s method here : Matlab code help on Euler's Method - MATLAB Answers - MATLAB Central (mathworks.com)

Make sure you set all the parameters as needed for your case. I would recommend running the above code on MATLAB R2023a for the best results.

Hope this helps!