I'm trying to solve this system of ODE's describing a mechanical spring model.
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (1)
Alan Stevens
am 2 Okt. 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 Kommentare
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!