工作日志excel php程序,科学网—matlab工作日志 - 钱磊的博文

20140916

Hello World!

disp('Hello World!')

20140324

Matlab处理fits文件一直是一个老大难问题，下面的连接提供了一套解决方案。

例子：

data=fits_read('myFITS.fits');header=fitsinfo('myFITS.fits');fits_write('myFITS_copy.fits', data, header);

20130906

ode45的精度可以通过odeset控制，具体如下。

ode = @(t,y) search3body(t,y);

options=odeset('reltol',1e-15,'abstol',1e-15)

[t,y] = ode45(ode, tspan, y0);

20130830

ode45的精度不知如何控制。步长可以调整。

tspan = 0:0.0001:10;

[t,y] = ode45(ode, tspan, y0);

20130827

运行了计算轨道的程序，熟悉了一下画图。

===========================================

% in search3body.m

function dydt = search3body(t,y)

% search3body Defines the equation of motion of three bodies of equal mass

% Copyright NAOC Qian Lei

% $Revision: 1.0 $  $Date: 2013/08/15 15:43:12 $

x1 = y(1);

y1 = y(2);

x2 = y(3);

y2 = y(4);

x3 = y(5);

y3 = y(6);

x1p = y(7);

y1p = y(8);

x2p = y(9);

y2p = y(10);

x3p = y(11);

y3p = y(12);

fx21 = (x2-x1)./(((x1-x2)^2+(y1-y2)^2)^(3/2));

fx31 = (x3-x1)./(((x1-x3)^2+(y1-y3)^2)^(3/2));

fy21 = (y2-y1)./(((x1-x2)^2+(y1-y2)^2)^(3/2));

fy31 = (y3-y1)./(((x1-x3)^2+(y1-y3)^2)^(3/2));

fx32 = (x3-x2)./(((x2-x3)^2+(y2-y3)^2)^(3/2));

fx12 = (x1-x2)./(((x2-x1)^2+(y2-y1)^2)^(3/2));

fy32 = (y3-y2)./(((x2-x3)^2+(y2-y3)^2)^(3/2));

fy12 = (y1-y2)./(((x2-x1)^2+(y2-y1)^2)^(3/2));

fx13 = (x1-x3)./(((x3-x1)^2+(y3-y1)^2)^(3/2));

fx23 = (x2-x3)./(((x3-x2)^2+(y3-y2)^2)^(3/2));

fy13 = (y1-y3)./(((x3-x1)^2+(y3-y1)^2)^(3/2));

fy23 = (y2-y3)./(((x3-x2)^2+(y3-y2)^2)^(3/2));

dydt = [x1p;...

y1p;...

x2p;...

y2p;...

x3p;...

y3p;...

fx21+fx31;...

fy21+fy31;...

fx32+fx12;...

fy32+fy12;...

fx13+fx23;...

fy13+fy23;...

];

===========================================

===========================================

% in test search3body.mtspan = [0, 100];

x1 = -1;

x2 = 1;

x3 = 0;

y1 = 0;

y2 = 0;

y3 = 0;

x1p = 0.1;

y1p = 0.2;

x2p = x1p;

x3p = -2*x1p;

y2p = y1p;

y3p = -2*y1p;

y0 = [x1; y1; x2; y2; x3; y3; x1p; y1p; x2p; y2p; x3p; y3p];

ode = @(t,y) search3body(t,y);

[t,y] = ode45(ode, tspan, y0);

% Plot of the solution

subplot(2,1,1)

hold on;

plot(y(:,1),y(:,2),'Color','red');

plot(y(:,3),y(:,4),'Color','green');

plot(y(:,5),y(:,6),'Color','blue');

xlabel('x')

ylabel('y')

title('3 body trajectories')

hold off;

format long e

sta = 1;

rpf = -log10(sqrt((y(sta:end,1)-x1).^2+(y(sta:end,2)-y1).^2+(y(sta:end,3)-x2).^2+(y(sta:end,4)-y2).^2 ...

+(y(sta:end,5)-x3).^2+(y(sta:end,6)-y3).^2+(y(sta:end,7)-x1p).^2+(y(sta:end,8)-y1p).^2 ...

+(y(sta:end,9)-x2p).^2+(y(sta:end,10)-y2p).^2 ...

+(y(sta:end,11)-x3p).^2+(y(sta:end,12)-y3p).^2));

subplot(2,1,2)

plot(t,rpf);

xlabel('t')

ylabel('r.p.f.')

title('3 body return proximity function')

===========================================

20130814

使用官网上的示例程序重新熟悉了一下用matlab求解微分方程。将微分方程写成一个函数保存在同名m文件中，然后用ode45进行求解。

============================================