UL/FRI/VSP-RI/NUM/Vaje/2006-11-23:B-D

računanje orbite kometov

Besedilo naloge na:

• http://matematika.fri.uni-lj.si/Numeri%C4%8DneMetode/Vaje/GibanjeKometa

Keplerjevi zakoni zakoni po katerih se komet giblje:

Nebesna telesa se gibljejo po elipsi. potrebno bo opisati elipso v parametrični obliki

Dobljeni podatki za elipso

• polos a
• ekscentričnost ε

ekscentričnost pove kakšno je razmerje polosi

Potrebujemo enačbo elipse, da bomo lahko narisali tir kometa.

 fi = linspace(0,2*pi);
> epsilon=0.976
epsilon = 0.97600
> a=17.788
a = 17.788
> r=sqrt(a^2*(1-epsilon^2)./(1 - epsilon^2 *cos(fi).^2));

> plot(fi,r)

(dobimo UU)

> x=r.*cos(fi);
> y=r.*sin(fi);

> plot(x,y)

(dobimo elipso)

Dorišemo gorišče elipse

> hold on > plot (a*epsilon,0,'x')

Newtonova metoda

• sestavimo zaporedje, ki konvergira proti rešitvi

zaporedje pri izbranem pravem začetnem približku konvergira proti nič

> f=inline('x-cos(x)','x')
f =

f(x) = x-cos(x)

> f(1)
ans = 0.45970
> f(pi)
ans = 4.1416
> df=inline('1+sin(1)','x')
df =

f(x) = 1+sin(1)

Konvergenca:

> x=3
x = 3
> x=x-f(x)/df(x)
x = 0.833258015185787
> x=x-f(x)/df(x)
x = 0.745941937122599
> x=x-f(x)/df(x)
x = 0.739700748605542
> x=x-f(x)/df(x)
x = 0.739141173485928
> x=x-f(x)/df(x)
x = 0.739090240925718
> x=x-f(x)/df(x)
x = 0.739085598802360
> x=x-f(x)/df(x)
x = 0.739085175655634
> x=x-f(x)/df(x)
x = 0.739085137083814
> x=x-f(x)/df(x)
x = 0.739085133567807
> x=x-f(x)/df(x)
x = 0.739085133247306
> x=x-f(x)/df(x)

Ko ponavljamo funkcijo z začetnim približkom 3 funkcija skonvergira proti 0.73908...

Funkcija newton.m

Toleranca nam pove na koliko decimalk natančno želimo rešitev.

function [x,it]=newton(f,df,x0,param,tol,maxit)
%NEWTON resuje nelinearno enacbo f(x)=0 z Newtonovo metodo
%[x,it]=newton(f,df,x0,tol,maxit)
%x je priblizek za resitev, it stevilo korakov
%f je funkcija, df odvod, x0 zacetni priblizek, tol toleranca
%maxit maksimalno dovoljeno stevilo operacij
%param je vektor dodatnih parametrov za funkciji f in df

%dokler je napaka vecja od tolerance in se nismo dosegli max iteracij ponavljamo
it=0;
napaka=tol-1;
x=x0;

while (abs(napaka)>tol)&(it<maxit)
	it=it+1;	
	popravek=x-feval(f,x,param)/feval(df,x,param);
	x=x-popravek;
	napaka=abs(popravek);
end

feval

> sin(1)
ans = 0.841470984807897
> feval('sin',1)
ans = 0.841470984807897

uporaba

> [x,it]=newton(f,df,3,[],1e-3,100)
x = -0.81501
it = 10

Reševanje enačb:

x³ + 2x² − 1 = 0

> f = inline('x^3+2*x^2-1','x','param')
f =

f(x, param) = x^3+2*x^2-1

> df = inline('3*x^2+4*x','x','param')
df =

f(x, param) = 3*x^2+4*x

>[x,it]=newton(f,df,3,[],1e-3,100)
x = -0.22673
it = 100

Funkcija komet.m

function [fi,r]=komet(t,T,e,a)
% funkcija [fi,r]=komet(t,T,e,a) izracuna kot in oddaljenost kometa v trenutku t
% T je perioda kometa
% e ekscentricnost tira
% a vecja polos tira

%kot
% resujemo enacbo fi-e^2*sin(fi)-2*pi*t/T=0
% param=[e,t,T]
f=inline('fi-param(1)^2*sin(fi)-2*pi*param(2)/param(3)=0','fi','param');
f=inline('1-param(1)^2*cos(fi)','fi','param');
% resimo enacbo
fi = newton(f,df,2*pi*t/T,[e,t,T],le-10,100);
%oddaljenost
r = sqrt(a^2*(1-e^2)/(1-e^2*cos(fi)^2))