(Take your time and answer no dumb solutions or chat gpt codes, please. Please run the code before submitting)

Use the same variable names and write the 3 functions F = Force(x,ks,kc,l0,S) ; function l0 = Compute_Rest_Length(x) ; function z = Equations_Of_Motion(xv,m,ks,kc,l0,S) ; in octave. And answer the 6 questions below. write only the three correct codes asked above in 3 different steps. Then solve the 6 problems in different steps. Please folle the sections as well. Thanks!

#!/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=j. % S = matrix indicating which springs exist % S(i,j) = 1 if there is a spring between particles i and j (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 (i0 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

% 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