Please use the same variables and only write the TODO part

#!/usr/bin/env octave

% Only modify this file where you see a "% TODO"

% Functions that do not contain this should not be modified.

% If you modify them for debugging, please remove/comment those modifications.

% This line tells octave that this is not a function file.

1;

% Workaround fltk bug. You might need to comment this line out.

graphics_toolkit("gnuplot");

% Meanings of commonly encountered variables:

% n = number of particles

% x = 2*n dimensional column vector with the positions of the particles

% v = 2*n dimensional column vector with the velocities of the particles

% Note that the ordering is always x1 y1 x2 y2 x3 y3 ...

% ks = spring constant

% kc = penalty collision strength

% l0 = matrix of spring rest lengths.

% l0(i,j) = initial length of spring between particles i and j (i

% Note that l0(i,j) should not be used if i>=j.

% S = matrix indicating which springs exist

% S(i,j) = 1 if there is a spring between particles i and j (i

% Note that S(i,j) should not be used if i>=j.

% m = mass of the particles (all particles have the same mass)

% Computes the total potential energy for the system.

function E = Potential_Energy(x,ks,kc,l0,S)

n = rows(S);

E = 0;

% For each pair of particles (i

for i = 1:n

xi = x(2*i-1:2*i,1); % Position of particle i

for j = i+1:n

% If there is a spring between particles i and j

if S(i,j)>0

xj = x(2*j-1:2*j,1); % Position of particle j

% Compute the potential energy of the spring.

E = E + .5*ks/l0(i,j)*(norm(xi-xj)-l0(i,j))^2;

end

end

% Compute the potential energy of the penalty collision force for particle i

% If the particle is inside the radius 2 circle centered at the origin,

% then there is no collision (no force, zero potential energy).

% Otherwise, the energy rises quadratically with distance from the circle.

E = E + .5*kc*max(norm(xi)-2,0)^2;

end

end

% Computes the total energy for the system.

function E = Total_Energy(x,v,m,ks,kc,l0,S)

n = rows(S);

KE = 0;

% Kinetic energy for particle = 1/2 m ||v||^2

for i = 1:n

vi = v(2*i-1:2*i,1);

KE = KE + .5*m*(vi'*vi);

end

% Total energy is kinetic + potential

E = KE + Potential_Energy(x,ks,kc,l0,S);

end

% Computes the total force for all of the particles in the system. x is a 2*n

% dimensional vector. On exit, F should be a 2*n dimensional vector containing

% the total force on each of the n particles. The force can be deduced from the

% potential energy. In particular, F(k) is the negative partial derivative of PE

% with respect to x(k). Here k=1..2*n and PE is the quantity computed by

% Potential_Energy. Treat the Potential_Energy function as a regular math

% function, which is a function of its input x.

function F = Force(x,ks,kc,l0,S)

n = rows(S);

F = zeros(2*n,1);

% TODO

end

% Computes the total momentum p and angular momentum L for the system

function [p L] = Momentum(x,v,m)

n = rows(x)/2;

p = zeros(2,1); % 2D vector; total momentum

L = 0; % Scalar; total angular momentum

for i = 1:n

xi = x(2*i-1:2*i,1);

vi = v(2*i-1:2*i,1);

L = L + det([xi m*vi]);

p = p + m*vi;

end

end

% This function should compute the lengths of the springs. Don't worry about

% whether a spring actually exists; just compute the pairwise distances for all

% pairs. You only need to fill l0(i,j) if i

% be used. l0 should be an (n x n) matrix.

function l0 = Compute_Rest_Length(x)

n = rows(x)/2;

l0 = zeros(n);

% TODO

end

% These are the equations of motion. We are trying to solve the differential

% equations:

% x' = dx;

% v' = dv;

% Here, you need to figure out what dx and dv should be. Note that this is not

% a discretization. We have no notion of x^(n+1) or x^n. We are not trying to

% do forward Euler or backward Euler. We are just trying to figure what the

% differential equation is that we need to solve. The equation will be solved

% for you. You may use the routines above to help you.

function z = Equations_Of_Motion(xv,m,ks,kc,l0,S)

n = rows(S);

x = xv(1:2*n);

v = xv(2*n+1:4*n);

dx = zeros(2*n,1);

dv = zeros(2*n,1);

% TODO

z = [dx ; dv];

end

% This function runs the simulations and plots the results.

% x0, v0 are 2*n dimensional arrays with initial positions and velocities for the particles.

% T is the final time.

% N is the number of simulation steps. There will be N+1 (initial + N steps) frames generated.

% m is the particle mass.

% ks is the spring stiffness.

% kc is the penalty collision stiffness.

% S indicates which springs exist (see comment above)

% label is a string to give output files different names.

function Analyze(x0,v0,T,N,m,ks,kc,S,label)

n = rows(S);

% Compute the spring rest lengths.

l0 = Compute_Rest_Length(x0);

% We will plot the position of the particles at times ts = 0, T/N, 2*T/N, ..., T.

ts = (0:N)'/N*T;

% Let octave solve the differential equation for us. This routine is actually

% extremely smart. It figures out what time step sizes to use, how many steps

% to take, what discretization to use, etc. It ensures that the simulation is

% stable and accurate for us. Note that although we only ask for N steps, it

% will often compute many more than this in order to achieve its accuracy

% requirements. XV will be a (N+1)x(4*n) matrix. The first index is the time

% step (ts has N+1 entries). The second index is the degree of freedom. The

% first 2*n entries are positions, the last 2*n entries are velocities.

tic

XV = lsode(@(xv,t) Equations_Of_Motion(xv,m,ks,kc,l0,S), [x0 ; v0], ts);

toc

% Get just the positions and clamp them to the plotting area.

X=max(min(XV(:,1:2*n),3),-3);

% Plot the collision area, which is a circle centered at the origin with radius 2.

o=(0:100)/100*2*pi;

plot(2*cos(o),2*sin(o),'-k');

axis([-3,3,-3,3]);

axis('square');

% Draw all the springs at all times, using a different color for each time.

C=jet(N+1);

hold on

for t = 0:N

for i = 1:n

xi = X(t+1,2*i-1:2*i);

for j = i+1:n

if S(i,j)>0

xj = X(t+1,2*j-1:2*j);

plot([xi(1);xj(1)],[xi(2);xj(2)],'-or','color',C(N-t+1,:));

end

end

end

end

% Dump the plot to disk.

hold off

print("-dpng","-color",['deformable-' label '.png'],'-r200');

pause(1);

% Plot the total energy, momentum, and angular momentum. The values are

% normalized to the range [-1,1] so that they can be plotted together in the

% same plot.

TE=zeros(N+1,1);

px=zeros(N+1,1);

py=zeros(N+1,1);

L=zeros(N+1,1);

for t = 0:N

x = XV(t+1,1:2*n)';

v = XV(t+1,2*n+1:4*n)';

TE(t+1) = Total_Energy(x,v,m,ks,kc,l0,S);

[pt,Lt] = Momentum(x,v,m);

px(t+1)=pt(1);

py(t+1)=pt(2);

L(t+1)=Lt;

end

TE = TE / max(TE);

norm_p=max(max(max(abs(px)),max(abs(py))),1e-15);

px = px / norm_p;

py = py / norm_p;

L = L / max(max(abs(L)),1e-15);

plot(ts,TE,'-or',ts,px,'-og',ts,py,'-ob',ts,L,'-om');

title(["Total energy (Test " label ")"]);

legend({"energy","momentum(x)","momentum (y)","angular momentum"},"location","southeast");

print("-deps","-color",['energy-' label '.eps']);

pause(1);

% For each time step, make a plot showing just the configuration at that

% time. Save each frame to a separate file.

for t = 0:N

plot(2*cos(o),2*sin(o),'-k');

axis([-3,3,-3,3]);

axis('square');

C=jet(N+1);

hold on

for i = 1:n

xi = XV(t+1,2*i-1:2*i);

for j = i+1:n

if S(i,j)>0

xj = XV(t+1,2*j-1:2*j);

plot([xi(1);xj(1)],[xi(2);xj(2)],'-or');

end

end

end

hold off

print("-dpng","-color",sprintf('frame-%04d.png',t+1),'-r200');

end

% Combine the per-frame simulation results into a video.

% Note that you may need to install ffmpeg for this line to work.

% This works on Linux and Mac, but it might not work on Windows.

% You can generate the movie from the frame-*.png files in another

% way if this does not work but you want to see the movie.

system(['ffmpeg -y -i frame-%04d.png -pix_fmt yuv420p -vcodec libx264 -b 4096k deformable-' label '.mov']);

end

% Number of time steps per second. Bigger numbers give you more output, which

% is rather like slow motion. Smaller numbers make things faster, but the video

% also goes by quicker.

N = 60;

% First test case. Triangle with vertices at (0,0) (1,0) and (0,1). Vertices

% are all moving with velocity (-1,-2), so the triangle is initially uniformly

% translating.

Analyze([0; 0; 1; 0; 0; 1],

[-1; -2; -1; -2; -1; -2],

2, 2*N, 10, 1e4, 1e5,

[0 1 1 ; 0 0 1 ; 0 0 0],

'A');

% Six triangles making a regular hexagon shape. This object starts off spinning

% and translating.

S = [

0 1 1 1 1 1 1 ;

0 0 1 0 0 0 1 ;

0 0 0 1 0 0 0 ;

0 0 0 0 1 0 0 ;

0 0 0 0 0 1 0 ;

0 0 0 0 0 0 1 ;

0 0 0 0 0 0 0

];

x0 = [0 cos((0:5)/6*2*pi); 0 sin((0:5)/6*2*pi)];

v0 = [0 sin((0:5)/6*2*pi); 0 -cos((0:5)/6*2*pi)]+2;

Analyze(reshape(x0,14,1), reshape(v0,14,1), 4, 4*N, 10, 1e4, 1e4, S, 'B');

% QUESTION 1: The plot of total energy does not change much. (If your energy

% plot does change significantly, you have a bug.) Why do you think this is?

% TODO

% QUESTION 2: Momentum (x and y) stays constant for a while, then rapidly

% changes value, then stays constant at that value for a while, etc. Why does

% it stay constant for a significant portion of the time? Why does its value

% sometimes change rapidly?

% TODO

% QUESTION 3 (Extra credit): Angular momentum stays constant for the entire

% simulation. Why? To get extra credit, your explanation must account for the

% following observations: (1) momentum changes sometimes but angular momentum

% never does and (2) if the center of the collision circle is moved, the angular

% momentum will no longer be constant.

% TODO

% QUESTION 4. Make the spring constant ks larger by a factor of 10. What

% qualitative change is observed in the results?

% TODO

% QUESTION 5. Make the collision constant kc smaller by a factor of 10. What

% qualitative change is observed in the results?

% TODO

% QUESTION 6: Make the mass smaller by a factor of 10. What qualitative change

% is observed in the results?

% TODO