Mengapa solusi numerik dari ODE menjauh dari keseimbangan yang tidak stabil?

15

Saya ingin mensimulasikan perilaku sistem seperti pendulum ganda. Sistem ini adalah manipulator robot 2 derajat kebebasan yang tidak digerakkan dan oleh karena itu, akan berperilaku seperti pendulum ganda yang dipengaruhi oleh gravitasi. Satu-satunya perbedaan utama dengan pendulum ganda adalah bahwa ia terdiri dari dua benda kaku dengan massa dan sifat inersia di pusat massa mereka.

Pada dasarnya, saya memprogram di ode45bawah Matlab untuk memecahkan sistem ODE dari jenis berikut:

[10000M110M1200100M120M22][x˙1x˙2x˙3x˙4]=[x2V1G1x4V2G2]

di mana x1 adalah sudut tubuh pertama sehubungan dengan horizontal, x2 adalah kecepatan sudut tubuh pertama; x3 adalah sudut tubuh kedua sehubungan dengan tubuh pertama, dan x4 adalah kecepatan sudut tubuh kedua. Semua koefisien ditentukan dalam kode berikut, dalam rhsdan fMassfungsi yang saya buat.

clear all
opts= odeset('Mass',@fMass,'MStateDependence','strong','MassSingular','no','OutputFcn',@odeplot);
sol = ode45(@(t,x) rhs(t,x),[0 5],[pi/2 0 0 0],opts);

function F=rhs(t,x)
    m=[1 1];
    l=0.5;
    a=[0.25 0.25];
    g=9.81;
    c1=cos(x(1));
    s2=sin(x(3));
    c12=cos(x(1)+x(3));
    n1=m(2)*a(2)*l;
    V1=-n1*s2*x(4)^2-2*n1*s2*x(2)*x(4);
    V2=n1*s2*x(2)^2;
    G1=m(1)*a(1)*g*c1+m(2)*g*(l*c1+a(2)*c12);
    G2=m(2)*g*a(2)*c12;

    F(1)=x(2);
    F(2)=-V1-G1;
    F(3)=x(4);
    F(4)=-V2-G2;
    F=F';     
end

function M=fMass(t,x)
    m=[1 1];
    l=0.5;
    Izz=[0.11 0.11];
    a=[0.25 0.25];
    c2=cos(x(3));
    n1=m(2)*a(2)*l;
    M11=m(1)*a(1)^2+Izz(1)+m(2)*(a(2)^2+l^2)+2*n1*c2+Izz(2);
    M12=m(2)*a(2)^2+n1*c2+Izz(2);
    M22=m(2)*a(2)^2+Izz(2);
    M=[1 0 0 0;0 M11 0 M12;0 0 1 0;0 M12 0 M22];
end

Perhatikan bagaimana saya mengatur kondisi awal x1 (sudut tubuh pertama sehubungan dengan horizontal) sehingga sistem mulai dalam posisi yang sepenuhnya vertikal. Dengan cara ini, karena hanya gravitasi yang bertindak, hasil yang jelas adalah bahwa sistem tidak boleh bergerak sama sekali dari posisi itu.

CATATAN: dalam semua grafik di bawah ini, saya merencanakan solusi x1 dan x3 sehubungan dengan waktu.

ODE45

Ketika saya menjalankan simulasi selama 6 detik ode45, saya mendapatkan solusi yang diharapkan tanpa masalah sama sekali, sistem tetap berada di tempatnya dan tidak bergerak:

masukkan deskripsi gambar di sini

Namun, ketika saya menjalankan simulasi selama 10 detik, sistem mulai bergerak secara tidak masuk akal:

masukkan deskripsi gambar di sini

ODE23

Saya kemudian menjalankan simulasi ode23untuk melihat apakah masalah tetap ada. Saya berakhir dengan perilaku yang sama, hanya kali ini perbedaan dimulai 1 detik kemudian:

masukkan deskripsi gambar di sini

ODE15s

Saya kemudian menjalankan simulasi dengan ode15suntuk melihat apakah masalah tetap ada dan tidak, sistem tampaknya stabil bahkan selama 100 detik:

masukkan deskripsi gambar di sini

Kemudian lagi, ode15shanya urutan pertama dan perhatikan bahwa hanya ada beberapa langkah mengintegrasikan. Jadi saya menjalankan simulasi lain dengan ode15sselama 10 detik tetapi MaxStepukuran0.01 untuk meningkatkan presisi, dan sayangnya, ini mengarah ke hasil yang sama dengan keduanyaode45danode23.

masukkan deskripsi gambar di sini

Biasanya, hasil nyata dari simulasi ini adalah bahwa sistem tetap pada posisi awalnya karena tidak ada yang mengganggu itu. Mengapa perbedaan ini terjadi? Apakah itu ada hubungannya dengan fakta bahwa jenis sistem ini pada dasarnya kacau? Apakah ini perilaku normal untuk odefungsi di Matlab?

jrojasqu
sumber
Selain persamaan, saya pikir skema juga akan banyak membantu untuk memahami pertanyaan.
nicoguaro
Jika menurut Anda itu sesuai, Anda dapat menerima salah satu jawaban (ada tombol hijau).
Ertxiem - mengembalikan Monica
Anda tidak mengatakannya, tetapi Anda tampaknya merencanakan x1dan x3. (Masukkan komentar kering tentang grafik tanpa legenda atau deskripsi.) Cobalah merencanakan logaritma (nilai absolut dari) x2dan x4.
Eric Towers

Jawaban:

15

Saya pikir dua poin utama telah dibuat oleh Brian dan Ertxiem: nilai awal Anda adalah keseimbangan yang tidak stabil dan fakta bahwa perhitungan numerik Anda tidak pernah benar-benar tepat memberikan gangguan kecil yang akan membuat ketidakstabilan menendang.

Untuk memberikan sedikit lebih detail bagaimana ini terjadi, pertimbangkan masalah Anda dalam bentuk masalah nilai awal umum

y˙(t)=M1f(t,y(t))
y(t)=(x1(t),x2(t),x3(t),x4(t))

f(t,y(t))=[x2V1G1x4V2G2]

f(0,y0)=0y˙(0)=0f~

Dalam kode Anda, Anda dapat mengujinya dengan komputasi

norm(rhs(0,[pi/2 0 0 0]))

yang memberikan 6.191e-16 - jadi hampir tetapi tidak sepenuhnya nol. Bagaimana hal itu berdampak pada dinamika sistem Anda?

fy0y~0 .

Selain itu, dalam waktu yang sangat singkat, solusi untuk sistem Anda terlihat seperti solusi dari sistem linierisasi

y˙(t)=f(0,y0)+f(0,y0)(y(t)y0)=f(0,y0)(y(t)y0)

ffrhsy0d(t):=y(t)y0d

d˙(t)=f(0,y0)d(t).

Saya tidak dapat diganggu untuk menghitung Jacobian dengan tangan jadi saya menggunakan diferensiasi otomatis untuk mendapatkan perkiraan yang baik:

J:=f(0,y0)=[01009.8102.4525000012.452502.45250]

so that your equation becomes

d˙(t)=Jd(t),d(0)=y~0y0

Now we need one final step: we can compute an eigenvalue decomposition of the Jacobian such that

J=QDQ1

where D is a diagonal matrix containing the eigenvalues of J and Q are orthogonal matrices, representing coordinate transformations. We can then transform the equation for d into an equation for e(t):=Q1d(t) which reads

e˙(t)=De(t),e(0)=Q1d0.

Because D is diagonal, these are effectively four independent equations

e˙i(t)=λiei(t),ei(0)=ith component of Q1d0

with i=1,2,3,4. If you compute the eigenvalues, you'll find that the largest one is λ1=3.2485. Therefore,

e1(t)=e1(0)e3.2485t.

Now if arithmetic in your computer were exact, you'd have d(0)=0, thus e(0)=Q1d(0)=0 and thus e1(0)=0 and nothing would happen. But since this is not the case, you have a small but finite e1(0) which gets exponentially amplified. Hence the rapid deviation from the equilibrium in your solution.

Daniel
sumber
16

Note that π/2 is represented in double precision format in a way that is not exactly equal to π/2. It's only accurate to about 15 digits. Thus you're starting every so slightly away from the equilibrium position. Since the equilibrium is unstable, it will eventually start moving.

Brian Borchers
sumber
4
If you monitor the state variables carefully (by looking at the values printed out in scientific notation), you should be able to see the initial very slow movement away from equilibrium.
Brian Borchers
This makes sense and indeed, when I start the system in a downwards vertical position (being a point of stable equilibrium), the system does not move at all, at least for a simulation of 1000 seconds which I consider a very long period of time.
jrojasqu
6
If you want to solve the problem, change your parametrization so that x1 is the angle with respect to the vertical axis. Then in your equations you'll need to compute sin(0) and cos(0) which are exact, rather than sin(pi/2) and cos(pi/2). In this way I think you can guarantee that rhs(t,[0,0,0 0] == [0,0,0,0], which means your system won't move at all.
Federico Poloni
@jrojasqu just to make explicitly clear what has been said: your variables are represented by double precision floats, which can only represent a finite set of numbers, and π/2 is not one of them for several reasons, the simplest being that it's irrational and so has an infinite number of digits, but floats have a finite number of digits. Zero is a member of this set and so it can be exactly represented.
llama
1
Note that if the straight down direction is not represented by θ=0, you should also see some movement where you would mathematically expect none, but in this case it will just be oscillating about the stable equilibrium with an amplitude on the order of the float precision ( 1016 for doubles)
llama
4

Look at the components of the forces calculated in your functions.

You will probably find they are never exactly zero, because as other answers have said, you can't represent the value of π exactly in computer arithmetic, and the routines that calculate trig functions are not exact either.

Eventually, the tiny forces (probably of order 1016 at the start) will move the system away from its unstable equilibrium position.

While the displacement of the system is still very small, all the calculations will lose a lot of precision through rounding errors (you are doing calculations similar to a=1.0; a=a+1016) so the amount of time before the system "topples over" in the simulation will depend on exactly what integration method you used, what time steps you requested, etc.

alephzero
sumber
4

The initial assumption was that the initial position was at a stable equilibrium (i.e., a minimum of the potential energy) with zero kinetic energy and the system started moving away from the equilibrium.
Since physically it can't happen (if we consider classical mechanics), two things came to my mind:

  1. The first one is that maybe the initial position is: both pendulums pointing upwards (π/2 instead of π/2?), which is a point of unstable equilibrium;

  2. The second one is that perhaps there is something wrong with the equations of movement (perhaps a typo somewhere?). Can you please write the equations explicitly? Perhaps you could plot the angular acceleration as a function of the initial position of each pendulum, assuming zero angular velocity to check if there is something weird.

Ertxiem - reinstate Monica
sumber
1
Indeed, I started the system in an upwards vertical position. Therefore, it is a point of unstable equilibrium. Brian Borcher's comment completes your answer by explaining the issue with π approximation which leads to the system eventually moving from that position.
jrojasqu
2
By the way, just for fun, if you wanted to keep the system at the unstable vertical position, you could change your origin of coordinates to have the angle equal to zero pointing upwards.
Ertxiem - reinstate Monica
@Ertxiem another option is to introduce small friction in the pins that would eat numerical errors.
svavil
@Ertxiem For the sake of fun, I tried changing the coordinate system so that the zero angle makes the system point upwards. It is really the best parameterization here. Obviously the system stays in the upwards position indefinitely. However oscillations still arise (minimal for 1000 seconds of simulation but there) at the stable equilibrium position (straight downwards) because then, there is a sin(π) to be calculated in the force deriving from the potential energy. So I insist in the fact that if I simulate this long enough, the system will start deviating also from that position.
jrojasqu
Since physically it can't happen – Given the insight that we are at an unstable equilibrium, I somewhat challenge this. Physical systems (without too much friction) do not stay in unstable equilibria. More generally, if you simulate real systems, you would want to avoid that it accidentally gets stuck in an unstable equilibrium (however it got there) – it’s a feature, not a bug. (There are some rare exceptions to this, such as the uninfected state in immunology, which is an unstable equilibrium that can be maintained.)
Wrzlprmft
0

You should search more about double pendulums: they are what we call "chaotic systems". Even though they behave following simple rules, starting from slightly different initial conditions, solutions diverge quite fast. Doing numerical simulations for this kind of systems is not easy. Take a look at the following video to get more insight into the problem.

For the simple or double pendulum you could write a formula for the total energy of the system. Assuming that friction forces are neglected, this total energy is preserved by the analytic system. Numerically this is a whole other issue.

Before trying the double pendulum, try the simple pendulum. You will notice that for Runge-Kutta methods of small order the energy of the system will grow in the numerical simulations, instead of remaining constant (this is what happens in your simulations: you get movement out of nothing). In order to prevent this, higher order RK methods could be used (ode45 is of order 4; RK of order 8 would work better). There are other methods called "symplectic methods" which are designed such that the numerical simulations conserve the energy. In general you should stop the simulation as soon as the energy increases significantly compared to your initialization.

And try to understand the simple pendulum before going to the double one.

Beni Bogosel
sumber
2
This is not a matter of the system being chaotic, though. You can have an unstable equilibrium in non-chaotic systems as well, e.g. single pendulum being "on its head", and it will exhibit the same behaviour described in the question.
Daniel
1
It is also not true that the energy increases for RKM of small order: implicit Euler ist first order and shows exactly the opposite behaviour.
Daniel
@BeniBogosel You mention symplectic methods which caught my attention because obviously, in my example, the energy is not conserved. However, could you indicate a specific symplectic method that could be implemented here?
jrojasqu
@jrojasqu why do you say that energy is not conserved on your system?
Ertxiem - reinstate Monica
@Ertxiem When I calculate the total Mechanical Energy of the system (Kinetic+Potential energies) with the outputs provided by ode45, I get a value that starts with zero, then grows with time. The value is very very small at the first seconds, but still, it consistently grows away from zero. I believe it is because of the issues that were addressed in the answers above (approximation of π, etc.).
jrojasqu