Plotting the Muller-Brown Potential and Trajectory Along It

Hey everyone,

I've been trying to plot the Muller-Brown potential using a contour plot and then simulation a trajectory along it. Here is the entirety of my code currently.

%known constants
A = [-200, -100, -170, 15];
a = [-1, -1, -6.5, 0.7];
b = [0, 0, 11, 0.6];
c = [-10, -10, -6.5, 0.7];
x0 = [1, 0, -0.5, -1];
y0 = [0, 0.5, 1.5, 1];

%commands to recreate the mapping of the PES
xx = linspace(-1.5, 1, 100);
yy = linspace(-0.5, 2, 100);
[X,Y] = meshgrid(xx,yy);


%Muller-Brown Potential Equation
W = A(1).*exp(a(1).*((X-x0(1)).^2) + b(1).*(X-x0(1)).*(Y-y0(1)) + c(1).*((Y-y0(1)).^2)) + A(2).*exp(a(2).*((X-x0(2)).^2) + b(2).*(X-x0(2)).*(Y-y0(2)) + c(2).*((Y-y0(2)).^2)) + A(3).*exp(a(3).*((X-x0(3)).^2) + b(3).*(X-x0(3)).*(Y-y0(3)) + c(3).*((Y-y0(3)).^2)) + A(4).*exp(a(4).*((X-x0(4)).^2) + b(4).*(X-x0(4)).*(Y-y0(4)) + c(4).*((Y-y0(4)).^2));
Z = W - min(W, [], 'all');
Z(Z >= 180) = NaN;

%plotting commands for the original PES
colormap('default')
contourf(X, Y, Z, 29)
hold on

%gradient equation
[DX, DY] = gradient(Z);
quiver(X, Y, DX, DY)
plot(1,0,'.', 'MarkerSize', 20) %plots starting point of the trajectory
hold off

%beginning of forward euler method to map trajectory on PES


h = 2*10^-5; %time step size
t = 2*10^4; %length of simulation
%N = t/h; %number of steps to be taken per simulation
N = 100; %shorter number of runs for practice
XX = zeros(1,100); %second value is the value of N, must be entered as an integer
YY = zeros(1,100);
YY(1) = 0; %initial condition for Y
XX(1) = 1; %initial condition for X

%equations for adjusting trajectory in x and y dimensions
%included for reference purposes currently
%dx = -DX + noise*sqrt(beta);
%dy = -DY + noise*sqrt(beta);

for i = 1:(N-1)
    XX(i+1) = XX(i) + h*f(XX(i));
    YY(i+1) = YY(i) + h*g(YY(i));
end

figure
colormap('default')
contourf(X, Y, Z, 29)
hold on
plot(XX,YY, 'o', 'MarkerSize', 5S)
hold off

function dx = f(DX)
    noise = 1;
    beta = 40;
    dx = -DX + noise*sqrt(beta);
end

function dy = g(DY)
    noise = 1;
    beta = 40;
    dy = -DY + noise*sqrt(beta);
end

The issue I am having is in the implementation of the forward Euler method. For some reason the values of XX and YY increase constantly by only the time step and I am really struggling to figure out how to fix that.

Thanks in advance if anyone can get this.

Answer 1

[Kavya Vuriti]

The values “XX” and “YY” are not increasing constantly according to your code. There is a slight difference in the increments. Set the output format to long fixed-decimal format to view this difference.

format long 
diff1=XX(3)-XX(2); 
diff2=XX(2)-XX(1);

The values diff1 and diff2 are different.

[Walter Roberson]

Zachary Raboin comments to Kavya Vuriti :

this didnt help