Giải phương trình chuyển động hệ sao kép

binary_star

%binary star system - He sao doi
clear all; close all; clc;
M1=2;M2=3; %M=1, mass of star 1 and 2
f=inline('[u(2); -20*pi^2*u(1)/(u(1)^2+u(3)^2 )^(3/2); u(4);...
-20*pi^2*u(3)/(u(1)^2+u(3)^2 )^(3/2)]','t','u');
opt=odeset('reltol',1e-9,'abstol',1e-12);
[t, u]=ode45(f,[0:0.001: .5],[-4/3;0;0;-4.0 ], opt);
%  [t, u]=ode45(f,[0 .5],[-4/3;0;0;-4.0 ], opt);  %can chu y

%chuyen dong cua he co khoi luong rut gon
x=u(:,1); y=u(:,3);
vx=u(:,2); vy=u(:,4);
r=sqrt(x.^2+y.^2); v= sqrt(vx.^2+vy.^2);
Ek=v.^2/2;%(Dong nang): v^2/2
Eu=-20*pi^2./r; %(The nang)
E=Ek+Eu;%(Nang luong toan phan)
Lz=x.*vy-y.*vx;% Mo men dg lg/m:Lz/m=x*vy-y*vx
figure(1)
subplot(3,1,1); plot(t,Ek,t,Eu,t,E,t,Lz);
legend('Ek(t)','Eu(t)','E(t)','Lz(t)'); title('Nang luong cua he rut gon')
subplot(3,1,2); plot( t, x, t, y);
legend('x(t)','y(t)');title('vi tri phu thuoc vao thoi gian cua he rut gon');
subplot(3,1,3); plot(t,vx,t,vy);
legend('vx(t)','vy(t)');title('van toc phu thuoc vao thoi gian cua he rut gon');

%thien the 1
x1=-M2/(M1+M2)*u(:,1);y1=-M2/(M1+M2)*u(:,3);
vx1=-M2/(M1+M2)*u(:,2);vy1=-M2/(M1+M2)*u(:,4);
%Nang luong thien the 1:
Ek1=0.5*M1*(vx1.^2+vy1.^2);
Eu1=-4*pi^2*M1*M2./sqrt(x.^2+y.^2);
E1=Ek1+Eu1;
Lz1=M1*(x1.*vy1-y1.*vx1);
figure(2);
subplot(3,1,1); plot(t,Ek1,t,Eu1,t,E1,t,Lz1);
legend('Ek1(t)','Eu1(t)','E1(t)','Lz1(t)');   title('nang luong cua thien the 1');
subplot(3,1,2); plot(t,x1,t,y1);
legend('x1(t)','y1(t)');title('vi tri phu thuoc vao thoi gian cua thien the 1');
subplot(3,1,3); plot(t,vx1,t,vy1);
legend('vx1(t)','vy1(t)');title('van toc phu thuoc vao thoi gian cua thien the 1');

%thien the 2
x2=M1/(M1+M2)*u(:,1); y2=M1/(M1+M2)*u(:,3);
vx2=M1/(M1+M2)*u(:,2); vy2=M1/(M1+M2)*u(:,4);
Ek2=0.5*M2*(vx2.^2+vy2.^2);
Eu2=-4*pi^2*M1*M2./sqrt(x.^2+y.^2);
E2=Ek2+Eu2;
Lz2=M2*(x2.*vy2-y2.*vx2);
figure(3);
subplot(3,1,1); plot(t,Ek2,t,Eu2,t,E2,t,Lz2);
legend('Ek2(t)','Eu2(t)','E2(t)','Lz2(t)');  title('nang luong cua thien the 2');
subplot(3,1,2); plot(t,x2,t,y2);
legend('x2(t)','y2(t)');title('vi tri phu thuoc vao thoi gian cua thien the 2');
subplot(3,1,3); plot(t,vx2,t,vy2);
legend('vx2(t)','vy2(t)');title('van toc phu thuoc vao thoi gian cua thien the 2');

%kiem tra xem nang luong tong cong co bang E1+E2 khong (tru di 1 lan the
%nang)
figure(4);
aE=Ek1+Ek2+Eu1;
plot(t,aE); ylim([-180 -160])

%Quy dao cua he sao
figure(5);
plot(x1,y1,'b--',x2,y2,'r--',u(:,1),u(:,3));hold on
h1=plot(x1(1),y1(1),'bo','markersize',4);
h2=plot(x2(1),y2(1),'ro','markersize',6);
h3=plot(u(1,1),u(3,1),'k.');
% set(h1, 'erase', 'xor');
% set(h2, 'erase', 'xor');
% set(h3, 'erase', 'xor');
axis equal;
axis([-1.5 1 -0.4 0.4])
for k=2:length(x1)
set( h1, 'xdata',x1(k),'ydata',y1(k) );
set( h2, 'xdata',x2(k), 'ydata',y2(k) );
set( h3, 'xdata',u(k,1), 'ydata',u(k,3) );
pause(0.001);
end

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Create a free website or blog at WordPress.com.

Up ↑

mlcvGru

Random thoughts on Machine Learning, Deep Learning and (sometimes) Computer Vision

LazyT

Góc nhỏ của Quyền

phanlan

we get what we give, sống là cho đi

Codeaholicguy

software engineer, team lead at #kobiton, blogger at @codeaholicguy

Maths 4 Physics & more...

Blog Toán Cao Cấp (M4Ps)

Vatlyvietnam's Blog

Thế Giới Song Song

Darren Wilkinson's research blog

Statistics, computing, data science, Bayes, stochastic modelling, systems biology and bioinformatics

Ông Xuân Hồng

Chia sẻ kiến thức và thông tin về Machine learning

Từ coder đến developer - Tôi đi code dạo

Lập trình viên giỏi không phải chỉ biết code

Computational Biology and Molecular Modelling

An interface between biology, chemistry and computer science

Moriator - I can do it!

Linux dễ dàng hơn bạn nghĩ!

VinaCode

Lập trình & Cuộc sống

Blog của Chiến

Học. Thực hành. Sáng tạo

Bespoke Blog

Science! Culture! Computational Engines!

Vuhavan's Blog

Just another WordPress.com weblog

%d bloggers like this: