2nd Order ODE Euler method not producing expected results
显示 更早的评论
Hi,
I am trying to solve a 2 degrees of freedom damped mass spring system using Euler's Method. I believe i am about 90% of the way there but the plot i am getting out is showing the displacement oscillations increasing exponentially(?) over time. I am pretty sure that the opposite should be happening.
Here is the starting matrix from which i've made my equations:
V and Z are 
respectively in order to solve with basic euler.
respectively in order to solve with basic euler.Here is my code:
close all;
clear;
clc
%Coefficients
m1=3;
m2=6;
k1=1;
k2=9;
k3=3;
c1=0.02;
c2=0.09;
c3=0.06;
% Step sizes
h=0.1;
t=0:h:100;
n=size(t);
N=n(2);
% Allocate array sizes by n steps
x1=zeros(n);
x2=zeros(n);
V=zeros(n);
Z=zeros(n);
% Initial Conditions
x1(1)=-7;
x2(1)=7;
V(1)=-3;
Z(1)=3;
% Create differential equations
dVdt=@(x1,x2,V,Z) 1/m1*(-(k1+k2)*x1+k2*x2-(c1+c2)*V+c2*Z);
dZdt=@(x1,x2,V,Z) 1/m2*(k2*x1-(k2+k3)*x2+c2*V-(c2+c3)*Z);
for i=1:N-1
V(i+1)=V(i)+h*dVdt(x1(i),x2(i),V(i),Z(i));
x1(i+1)=x1(i)+h*V(i);
Z(i+1)=Z(i)+h*dZdt(x1(i),x2(i),V(i),Z(i));
x2(i+1)=x2(i)+h*Z(i);
end
plot(t,x1, 'blue',t,x2,'red');
xlabel('t');
ylabel('x');
With the resulting plot ^^^
As you can see the x values become very large which definitely shouldn't be the case.
Can someone please point out the dumb thing I'm doing that's giving me this kind of result?
Thanks in advance!
采纳的回答
That is not uncommon with Euler. It can be an unstable method. At each step, a small amount of error gets added in. Eventually, the noise starts to dominate, if the step size is sufficiently large.
Better methods, (for example backwards Euler) were designed to avoid this problem.
One clue to this is to reduce the step size (h), and see if suddenly, the Euler method solution becomes well-behaved. For example, if I use h=0.0001 or 0.001, the result below is consistent. But you had h=0.1. This is apparently sufficiently large to push Euler into its unstable domain for this problem. Do some reading here:
m1=3;
m2=6;
k1=1;
k2=9;
k3=3;
c1=0.02;
c2=0.09;
c3=0.06;
% Step sizes
h=0.0001;
t=0:h:100;
n=size(t);
N=n(2);
% Allocate array sizes by n steps
x1=zeros(n);
x2=zeros(n);
V=zeros(n);
Z=zeros(n);
% Initial Conditions
x1(1)=-7;
x2(1)=7;
V(1)=-3;
Z(1)=3;
% Create differential equations
dVdt=@(x1,x2,V,Z) 1/m1*(-(k1+k2)*x1+k2*x2-(c1+c2)*V+c2*Z);
dZdt=@(x1,x2,V,Z) 1/m2*(k2*x1-(k2+k3)*x2+c2*V-(c2+c3)*Z);
for i=1:N-1
V(i+1)=V(i)+h*dVdt(x1(i),x2(i),V(i),Z(i));
x1(i+1)=x1(i)+h*V(i);
Z(i+1)=Z(i)+h*dZdt(x1(i),x2(i),V(i),Z(i));
x2(i+1)=x2(i)+h*Z(i);
end
plot(t,x1, 'blue',t,x2,'red');
xlabel('t');
ylabel('x');
My guess is this is probably what you were expected to trip over, why this problem was assigned to you. In general, Euler's method can be dangerous, certainly if you don't know what to watch out for.
2 个评论
James Webb
2023-3-18
John D'Errico
2023-3-18
编辑：John D'Errico
2023-3-18
Yep. They wanted you to have exactly this problem. I assume this would be a prelude to explaining exactly what we have told you, that forward Euler can be unstable when the step size is too large. The comparison to other solvers would show something strange, then making you endlessly search for what you did wrong. Thank your teacher. ;-)
This link does give a decent explanation of what happens, and why Euler can be unstable for some step sizes:
更多回答（1 个）
Not that bad.
%Coefficients
m1=3;
m2=6;
k1=1;
k2=9;
k3=3;
c1=0.02;
c2=0.09;
c3=0.06;
% Step sizes
h=0.00005;
t=0:h:100;
n=size(t);
N=n(2);
% Allocate array sizes by n steps
x1=zeros(n);
x2=zeros(n);
V=zeros(n);
Z=zeros(n);
% Initial Conditions
x1(1)=-7;
x2(1)=7;
V(1)=-3;
Z(1)=3;
% Create differential equations
dVdt=@(x1,x2,V,Z) 1/m1*(-(k1+k2)*x1+k2*x2-(c1+c2)*V+c2*Z);
dZdt=@(x1,x2,V,Z) 1/m2*(k2*x1-(k2+k3)*x2+c2*V-(c2+c3)*Z);
for i=1:N-1
V(i+1)=V(i)+h*dVdt(x1(i),x2(i),V(i),Z(i));
x1(i+1)=x1(i)+h*V(i);
Z(i+1)=Z(i)+h*dZdt(x1(i),x2(i),V(i),Z(i));
x2(i+1)=x2(i)+h*Z(i);
end
figure(1)
plot(t,x1, 'blue',t,x2,'red');
xlabel('t');
ylabel('x');
fun = @(t,y)[y(2);(-(k1+k2)*y(1)+k2*y(3)-(c1+c2)*y(2)+c2*y(4))/m1;y(4);(k2*y(1)-(k2+k3)*y(3)+c2*y(2)-(c2+c3)*y(4))/m2];
tspan = [0 100];
y0 = [-7 -3 7 3];
[T ,Y] = ode45(fun,t,y0,odeset('RelTol',1e-12,'AbsTol',1e-12));
figure(2)
plot(T,[Y(:,1),Y(:,3)])
figure(3)
plot(t,[Y(:,1).'-x1;Y(:,3).'-x2])
类别
在 帮助中心 和 File Exchange 中查找有关 Numeric Solvers 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
