Hi 현호 이
To extend the code for a system of two coupled oscillators to a system of three coupled oscillators, we need to adjust the system's dynamics and visualization. Changes include the addition of a third input for the initial position of the new oscillator and this will involve extending the arrays for position calculations (x3, y3) and adapting the spring visualization to include the third spring connecting the second and third masses.
Here is the adjusted code for your reference:
% Complete script for three coupled oscillators visualization
clf; clear; clc;
global x y z
a = input("Enter the initial value x0(-5<x0<5) : ");
b = input("Enter the initial value y0(-5<y0<5) : ");
c = input("Enter the initial value z0(-5<z0<5) : "); % Input for the third oscillator
osc(a, b, c); % Call to the osc function with three initial conditions
x1 = x + 5; y1 = 0;
x2 = y + 15; y2 = 0;
x3 = z + 25; y3 = 0; % Position adjustments for the third object
% Lines for springs visualization
t1 = linspace(0, x1(1), 100);
t2 = linspace(0, x2(1) - x1(1), 100);
t3 = linspace(0, x3(1) - x2(1), 100);
t4 = linspace(0, 30 - x3(1), 100);
bounds = [0 30 -5 5]; % Adjusted visualization bounds for third object
mysize = 70;
hold on
h = plot(x1(1), y1, "marker", ".", 'MarkerSize', mysize);
h2 = plot(x2(1), y2, "marker", ".", 'MarkerSize', mysize);
h3 = plot(x3(1), y3, "marker", ".", 'MarkerSize', mysize); % Marker for the third mass
% Spring visualizations
h4 = plot(t1, 1/10*sin(8*pi/x1(1)*t1), 'b');
h5 = plot(linspace(x1(1), x2(1), 100), 1/10*sin(8*pi/(x2(1)-x1(1))*t2));
h6 = plot(linspace(x2(1), x3(1), 100), 1/10*sin(8*pi/(x3(1)-x2(1))*t3), 'b'); % For the spring between mass 2 and 3
h7 = plot(linspace(x3(1), 30, 100), 1/10*sin(8*pi/(30-x3(1))*t4), 'b'); % For the right-most part beyond mass 3
axis(bounds); axis('equal');
pause(1);
for i_ = 1:length(x1)
grid on; axis(bounds); axis('equal');
set(h, 'XData', x1(i_))
set(h2, 'XData', x2(i_))
set(h3, 'XData', x3(i_)) % Updating position for the third mass
set(h4, 'XData', linspace(0, x1(i_), 100))
set(h4, 'YData', 1/10*sin(8*pi/x1(i_)*(linspace(0, x1(i_), 100))))
set(h5, 'XData', linspace(x1(i_), x2(i_), 100))
set(h5, 'YData', 1/10*sin(8*pi/(x2(i_)-x1(i_))*(linspace(0, x2(i_)-x1(i_), 100))))
set(h6, 'XData', linspace(x2(i_), x3(i_), 100))
set(h6, 'YData', 1/10*sin(8*pi/(x3(i_)-x2(i_))*(linspace(0, x3(i_)-x2(i_), 100))))
set(h7, 'XData', linspace(x3(i_), 30, 100))
set(h7, 'YData', 1/10*sin(8*pi/(30-x3(i_))*linspace(0, 30-x3(i_), 100)))
pause(0.03)
end
%---------------------------------------------------------------------------------------------------------------
function osc(a,b,c)
global x y z
x0=a; y0=b; z0=c;
tspan = linspace(0,15*pi,150);
m = 1;
k1 = 1; k2 = 3; k3 = 1; k4 = 1;
w1 = sqrt(k1/m); w2 = sqrt(k2/m); w3 = sqrt(k3/m); w4 = sqrt(k4/m);
A = [w1^2+w2^2, -w2^2, 0; -w2^2, w2^2+w3^2, -w3^2; 0, -w3^2, w3^2+w4^2];
[V,D] = eig(A);
Init = V * [x0; y0; z0];
X0 = Init(1); Y0 = Init(2); Z0 = Init(3);
Omega1 = sqrt(D(1,1)); Omega2 = sqrt(D(2,2)); Omega3 = sqrt(D(3,3));
X= X0* cos(Omega1*tspan);
Y= Y0* cos(Omega2*tspan);
Z= Z0* cos(Omega3*tspan);
x = X * V(1,1) + Y * V(2,1) + Z * V(3,1);
y = X * V(1,2) + Y * V(2,2) + Z * V(3,2);
z = X * V(1,3) + Y * V(2,3) + Z * V(3,3);
end
Hope this helps!
