Eulers for two degree spring system
显示 更早的评论
Hi there, I have to use basic Euler's to find and plot solutions to x1 and x2 over 0 < t < 30 where my step size is 0.001 with initial conditions: my conditions are x1(0) = 5, x2(0) = 1, v1(0) = 0,
Firstly I have multiplied the matrix attached out, subbed in v for x' and rearranged for x'' and renamed it as dv/dt where m = 1 c = 1 k = 2
i have obtained: dv/dt = 1 - (4(x1) - 2(x2)) - (2(v1) - (v2))
dv/dt = 1 - (4(x2) - 2(x1)) - (2(v2) - (v1))
I then tried to implement Euler's for the first equation on matlab t=zeros(30000,0); y=zeros(30000,0);
t(0)=0; x(0)=5; v(0)=0;
for i=1:30000
t(i+1)= t(i)+0.001;
x(i+1)= x(i)+0.001.*(1-(4*x(i)-2*x(i)) -(2*v(i)-v(i));
end
plot(t,x)
I'm getting an error in line 11, unbalances parenthesis and I'm not sure I'm really implementing the scheme correctly so any help would be greatly appreciated!
采纳的回答
James Tursa
2018-3-14
编辑：James Tursa
2018-3-14
You don't have enough equations in your code. This is a 4th order system (Two variables each as a 2nd order ODE gives 2x2 = 4 orders), so you will need 4 1st order equations. Using the nomenclature in the MATLAB doc for the ODE methods, let's pick y to be the state vector of your system (this will make it easier to compare your code with the ODE examples in the doc). Since you have a 4th order system, y will have 4 elements. Let's pick the elements to be the following:
y(1) = x1
y(2) = x2
y(3) = d(x1)/dt = x1dot = v1
y(4) = d(x2)/dt = x2dot = v2
The derivative of your state vector would then be
d(y(1))/dt = d(x1)/dt = y(3)
d(y(2))/dt = d(x2)/dt = y(4)
d(y(3))/dt = something
d(y(4))/dt = something
For those last two derivative equations, the "something" on the right hand side will be what you get by solving your two posted 2nd order equations for x1dotdot and x2dotdot and then substituting the appropriate y elements for the x1, x2, x1dot, and x2dot variables.
Then for your Euler iterations, at each step you will be propagating this 4 element vector. To do that, one way is to redefine our y vector as a matrix where each column is a y vector at a given time point. E.g.,
y(1,i) = x1 at time t(i)
y(2,i) = x2 at time t(i)
y(3,i) = x1dot at time t(i)
y(4,i) = x2dot at time t(i)
Then using this definition and the equations discussed above your MATLAB code would look something like this:
t = zeros(30001); % pre-allocate the time vector
y = zeros(4,30001); % pre-allocate the state vector
h = 0.001; % The stepsize
t(1) = 0; % Initial time (needs to start at index 1, not index 0)
y(1,1) = 5; % Initial x1 at time t(1) is contained in first column at y(1,1)
y(2,1) = 1; % Initial x2 at time t(1) is contained in first column at y(2,1)
y(3,1) = 0; % Initial v1 at time t(1) is contained in first column at y(3,1)
y(4,1) = 0; % Initial v2 at time t(1) is contained in first column at y(4,1)
for i=1:30000 % Note that the i=1 spot is already filled with initial conditions
t(i+1)= t(i) + h;
y(1,i+1) = y(1,i) + h * ( );
y(2,i+1) = y(2,i) + h * ( );
y(3,i+1) = y(3,i) + h * ( );
y(4,i+1) = y(4,i) + h * ( );
end
plot(t,y(1,:),t,y(2,:)); % y(1,:) is the x1 solution, and y(2,:) is the x2 solution
What you need to do is fill in the blanks inside the parentheses in the code above. Everything inside those parentheses will be an expression involving y(1,i), y(2,i), y(3,i), and y(4,i). The first two should be easy based on the above discussion. The last two will involve your solution to solving your two posted 2nd order equations. Give it a try, and if you get stuck go ahead and post how far you got and what you are stuck on.
3 个评论
gary stewart
2018-3-24
gary stewart
2018-3-24
James Tursa
2018-3-24
Tackle the first two equations first. From our my definition of y above and adding in the column index i we have this:
d(y(1,i))/dt = d(x1_i)/dt = y(3,i)
d(y(2,i))/dt = d(x2_i)/dt = y(4,i)
So that tells you what to put in the parentheses for your first two equations. You get this:
y(1,i+1) = y(1,i) + h * (y(3,i));
y(2,i+1) = y(2,i) + h * (y(4,i));
For the last two equations, you use the solution for the two 2nd order equations you did above, but again write it in terms of y(1,i), y(2,i), y(3,i) and y(4,i). E.g.,
d(y(3,i))/dt = d(x1dot at time t(i))/dt = x1dotdot at time t(i)
d(y(4,i))/dt = d(x2dot at time t(i))/dt = x2dotdot at time t(i)
Then the expressions for x1dotdot (which is v1dot) and x2dotdot (which is v2dot) is from your solution to the 2nd order equations. You just need to plug in y(1,i), y(2,i), y(3,i) and y(4,i) for the appropriate x1, x2, v1, v2 elements.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
