gistlib - create a script that plots the trajectory of the geodesic equation in matlab

To create a script that plots the trajectory of the geodesic equation in Matlab, we can follow these steps:

Define the geodesic equation that we want to solve. This may depend on the specific problem we are interested in, but a general form for the geodesic equation in curved spacetime is:


main.m
function d2x = geodesic_equation(t,x,gamma)
% x is a four-vector representing the position and velocity of the particle
% gamma is the metric tensor
% d2x is the second derivative of x wrt proper time tau

dtau = sqrt(gamma(x(1:3),x(1:3))*x(4:end)'*x(4:end)); % proper time interval
G = gamma(x(1:3),:);

d2x = [x(4:end);...
       -(1/2)*sum(sum(inv(G).*reshape(sum(reshape(permute(sum(bsxfun(@times,G,reshape(x(4:end),1,1,[]))),2),[],3),size(G,1),size(G,2)),3).*reshape(bsxfun(@times,G(:,:,find(size(G,3)==size(G,2))),x(4:end)'),[],1,size(G,2)),2),1)/dtau^2];
562 chars
11 lines

This function calculates the second derivative of a four-vector x wrt proper time tau, given the metric tensor gamma.

Choose the initial conditions for the particle. For example, we may set the initial position to [0;0;0] and the initial 4-velocity to [1;0;0;0]. We also need to specify the range of proper times over which we want to solve the geodesic equation.


main.m
x0 = [0;0;0;1;0;0;0]; % initial position and velocity
tau_span = [0 10]; % range of proper times
97 chars
3 lines

Choose the metric tensor for the spacetime we are interested in. For example, we may use the Schwarzschild metric for a non-rotating black hole:


main.m
M = 1; % mass of the black hole
gamma = @(r) [-(1-2*M./r) zeros(1,3);...
              zeros(3,1) eye(3)*(1-2*M./r)];
118 chars
4 lines

This defines a function gamma that takes a radial coordinate r and returns the 4x4 metric tensor evaluated at that point. Note that we have set units where c=G=1.

Solve the geodesic equation using Matlab's ODE solver ode45. This will give us the trajectory of the particle as a function of proper time.


main.m
options = odeset('RelTol',1e-12,'AbsTol',1e-15); % specify the ODE solver options
sol = ode45(@(t,x) geodesic_equation(t,x,gamma), tau_span, x0, options); % solve the ODE
171 chars
3 lines

Plot the trajectory of the particle using Matlab's plot3 function.


main.m
x = sol.y(1:3,:);
plot3(x(1,:),x(2,:),x(3,:));
xlabel('x'); ylabel('y'); zlabel('z');
86 chars
4 lines

Putting these steps together, we obtain the following Matlab script:


main.m
% Define the geodesic equation
function d2x = geodesic_equation(t,x,gamma)
dtau = sqrt(gamma(x(1:3),x(1:3))*x(4:end)'*x(4:end)); % proper time interval
G = gamma(x(1:3),:);
d2x = [x(4:end);...
       -(1/2)*sum(sum(inv(G).*reshape(sum(reshape(permute(sum(bsxfun(@times,G,reshape(x(4:end),1,1,[]))),2),[],3),size(G,1),size(G,2)),3).*reshape(bsxfun(@times,G(:,:,find(size(G,3)==size(G,2))),x(4:end)'),[],1,size(G,2)),2),1)/dtau^2];

% Set the initial conditions and range of proper time
x0 = [0;0;0;1;0;0;0];
tau_span = [0 10];

% Define the metric tensor
M = 1;
gamma = @(r) [-(1-2*M./r) zeros(1,3);...
              zeros(3,1) eye(3)*(1-2*M./r)];

% Solve the geodesic equation
options = odeset('RelTol',1e-12,'AbsTol',1e-15);
sol = ode45(@(t,x) geodesic_equation(t,x,gamma), tau_span, x0, options);

% Plot the trajectory
x = sol.y(1:3,:);
plot3(x(1,:),x(2,:),x(3,:));
xlabel('x'); ylabel('y'); zlabel('z');
909 chars
25 lines

This script defines the geodesic equation, sets the initial conditions and spacetime metric, solves the geodesic equation using ode45, and plots the resulting trajectory. Note that this is just one example of how to do this - the specific form of the geodesic equation will depend on the problem at hand, and the metric tensor may need to be adjusted accordingly.

create a script that plots the trajectory of the geodesic equation in matlab

similar matlab code snippets

related categories