I'm trying to solve this system of ODE's describing a mechanical spring model.

I made all the equations symbolic functions and am trying to use a for loop to use Eulers method to solve them, but im getting really large numbers.

clear all; clc;

%applied forces
P=[1100; 1800; 3300];

%spring constants
k=[4500; 1650; 1100; 2250; 550; 9300];

%friction coefficient 
b=50;

m=100;



syms f1(x1,x2,x3,x4) f2(x1,x2,x3,x4) f3(x1,x2,x3,x4) f4(x1,x2,x3,x4) t
sympref('FloatingPointOutput',true);
%PART 1
%use eulers method for "simulation" to solve odes
%for long time
%xn+1=xn+hfn
%h=step size aka time difference
%mass 1
f1(x1,x2,x3,x4)=P(1)-k(1)*x1-k(2)*(x1-x2)-k(4)*(x1-x3)-m*diff(x1,t,2)-b*diff(x1,t)==0;

%mass 2
f2(x1,x2,x3,x4)=P(2)+k(2)*(x1-x2)-k(3)*x2-k(6)*(x2-x4)-m*diff(x2,t,2)-b*diff(x2,t)==0;

%mass 3
f3(x1,x2,x3,x4)=P(3)+k(4)*(x1-x3)-k(5)*(x3-x4)-m*diff(x3,t,2)-b*diff(x3,t)==0;

%mass 4
f4(x1,x2,x3,x4)=k(6)*(x2-x4)+k(5)*(x3-x4)-m*diff(x4,t,2)==0;


%initial value x=0


step=0.001;
x1=zeros([1 100]);
x2=zeros([1 100]);
x3=zeros([1 100]);
x4=zeros([1 100]);

for i=2:100
    
x1(i)=x1(i-1)+f1(x1(i-1),x2(i-1),x3(i-1),x4(i-1))*step;
x2(i)=x2(i-1)+f2(x1(i-1),x2(i-1),x3(i-1),x4(i-1))*step;
x3(i)=x3(i-1)+f3(x1(i-1),x2(i-1),x3(i-1),x4(i-1))*step;
x4(i)=x4(i-1)+f4(x1(i-1),x2(i-1),x3(i-1),x4(i-1))*step;
end

Can anyone see where I'm going wrong?

Answer 1

[Alan Stevens]

Might be better to forget about symbolics, treat each 2nd order ode as two first order ode's and do the following:

%applied forces
P=[1100; 1800; 3300];

%spring constants
k=[4500; 1650; 1100; 2250; 550; 9300];

%friction coefficient 
b=50;

m=100;

% dx1dt = v1
% dv1dt = (P-k1*x1-k2*(x1-x2)-k4*(x1-x3)-b*v)/m

%mass 1
f1v = @(x1,x2,x3,x4,v1) (P(1)-k(1)*x1-k(2)*(x1-x2)-k(4)*(x1-x3)-b*v1)/m;   % dv1/dt
f1x = @(v1) v1;                                                            % dx1/dt

%mass 2
f2v = @(x1,x2,x3,x4,v2) (P(2)+k(2)*(x1-x2)-k(3)*x2-k(6)*(x2-x4)-b*v2)/m;
f2x = @(v2) v2;
%mass 3
f3v = @(x1,x2,x3,x4,v3)(P(3)+k(4)*(x1-x3)-k(5)*(x3-x4)-b*v3)/m;
f3x = @(v3) v3;

%mass 4
f4v = @(x1,x2,x3,x4)(k(6)*(x2-x4)+k(5)*(x3-x4))/m;
f4x = @(v4) v4;

%initial value x=0

step=0.001;
t = 0:step:20;
x1=zeros(1,numel(t)); v1 = x1;
x2 = x1; v2 = v1;
x3 = x1; v3 = v1;
x4 = x1; v4 = v1;

for i=2:numel(t)
x1(i) = x1(i-1) + f1x(v1(i-1))*step;    
v1(i)=v1(i-1)+f1v(x1(i-1),x2(i-1),x3(i-1),x4(i-1),v1(i-1))*step;
x2(i) = x2(i-1) + f2x(v1(i-1))*step;    
v2(i)=v2(i-1)+f2v(x1(i-1),x2(i-1),x3(i-1),x4(i-1),v2(i-1))*step;
x3(i) = x3(i-1) + f3x(v1(i-1))*step;    
v3(i)=v3(i-1)+f3v(x1(i-1),x2(i-1),x3(i-1),x4(i-1),v3(i-1))*step;
x4(i) = x4(i-1) + f4x(v4(i-1))*step;
v4(i) = v4(i-1) + f4v(x1(i-1),x2(i-1),x3(i-1),x4(i-1))*step;
end

subplot(2,2,1)
plot(t,x1),grid
xlabel('t'),ylabel('x1')
subplot(2,2,2)
plot(t,x2),grid
xlabel('t'),ylabel('x2')
subplot(2,2,3)
plot(t,x3),grid
xlabel('t'),ylabel('x3')
subplot(2,2,4)
plot(t,x4),grid
xlabel('t'),ylabel('x4')