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

조회 수: 2(최근 30일)
Caylyn MacDougall
Caylyn MacDougall 2021년 10월 2일
답변: Alan Stevens 2021년 10월 2일
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?

답변(1개)

Alan Stevens
Alan Stevens 2021년 10월 2일
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')

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by