How to make 3-mass-spring coupled oscillator system?

I've already made about two mass spring coupled oscillator system.

But, I;m trying to make about three, but it didn't work out. please help me to make code.

This code is for two objects

function osc(a,b)
global x y
x0=a ; y0=b;
tspan = linspace(0,15*pi,150) ;
m= 1 ;
k1= 1; k2=3 ; k3=1 ;
w1= sqrt(k1/m); w2= sqrt(k2/m) ; w3= sqrt(k3/m) ;
A= [w1^2+w2^2 , -w2^2 ; -w2^2, w2^2+w3^2] ;
[V,D] = eig(A) ;
Init = V * [x0 ;y0] ;
X0 = Init(1) ; Y0 = Init(2) ;
Omega1 = sqrt(D(1,1)); Omega2 = sqrt(D(2,2));
X= X0* cos(Omega1*tspan) ;
Y= Y0* cos(Omega2*tspan) ;
x = X * V(1,1) + Y* V(2,1) ;
y = X * V(1,2) + Y* V(2,2) ;
end
%---------------------------------------------------------------------------------------------------------------
clf; clear all; clc;
global x y
a=input("Enter the inital value x0(-5<x0<5) : ")
b=input("Enter the inital value y0(-5<y0<5) : ")
osc(a,b);
x1=x+5; y1=0; x2=y+15; y2=0;
t1=linspace(0,x1(1),100);
t2=linspace(0,x2(1)-x1(1),100);
t3=linspace(0,20-x2(1),100);
bounds = [0 20 -5 5];
mysize=70;
hold on
h = plot(x1(1),y1,"marker", ".", 'MarkerSize', mysize);
h2 = plot(x2(1),y1,"marker", ".", 'MarkerSize', mysize);
h3 = plot(t1,1/10*sin(8*pi/x1(1)*t1),'b');
h4 = plot(linspace(x1(1),x2(1),100),1/10*sin(8*pi/(x2(1)-x1(1))*t2));
h5 = plot(linspace(x2(1),20,100),1/10*sin(8*pi/(20-x2(1))*t3),'b');
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',linspace(0,x1(i_),100))
set(h3, 'YData',1/10*sin(8*pi/x1(i_)*(linspace(0,x1(i_),100))))
set(h4, 'XData',linspace(x1(i_),x2(i_),100))
set(h4, 'YData',1/10*sin(8*pi/(x2(i_)-x1(i_))*(linspace(0,x2(i_)-x1(i_),100))))
set(h5, 'XData',linspace(x2(i_),20,100))
set(h5, 'YData',1/10*sin(8*pi/(20-x2(i_))*linspace(0,20-x2(i_),100)))
pause(0.03)
end

Answer 1

[Vaibhav]

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!