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

조회 수: 2(최근 30일)
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 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