I'm trying to solve this system of ODE's describing a mechanical spring model.
2 views (last 30 days)
Show older comments
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?
0 Comments
Answers (1)
Alan Stevens
on 2 Oct 2021
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')
0 Comments
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!