Zadanie 1.3

near0=(10^(-10)); (*Ustalam początek NDSolve ≠ 0, by uniknąć błędu (1/0)*)
start = near0; (*Ustalam zakres obliczeń NDSolve*)
stop=20;
boundaryCond={θ[near0]==1,θ'[near0]==0}; (*Warunki brzegowe*)

(*Równanie Lane-Emdena dla różnych n*)
eqLE=θ''[ξ]+(2/ξ)θ'[ξ]+((θ[ξ])^n)==0;

Dla n = 0

n = 0 w przybliżeniu opisuje strukturę Ziemi (planety skalistej).

Clear[n,sol0,ξ0,ξ1];n=0;

sol0=eqLE;
sol0=NDSolve[{sol0,boundaryCond,WhenEvent[Re[θ[ξ]]==0,{ξ0=ξ,"StopIntegration"}]},θ,{ξ,start,stop}];ξ1=Re[θ'[ξ0]/. sol0[[1]]];
Plot[θ[ξ]/. sol0[[1]],{ξ,start,ξ0},PlotRange->All,AxesLabel->{"ξ","θ"},GridLines->Automatic]

Gęstość

ρ[ξ_]:=(Subscript[ρ, c]) ((θ[ξ])^n)/.sol0[[1]];
Plot[ρ[ξ]/.(Subscript[ρ, c]) ->1,{ξ,0,ξ0},PlotRange->All,AxesLabel->{"ξ","ρ /(Subscript[ρ, c])"},GridLines->Automatic]

Gęstość średnia:

⟨ρ⟩=((NIntegrate[ρ[ξ] 4 π ξ^2/.(Subscript[ρ, c])->1,{ξ,near0,ξ0}])/((4/3)π ((ξ0)^3)))(Subscript[ρ, c])

Kontrast gęstości:

((Subscript[ρ, c])/(⟨ρ⟩))

-((ξ0)/(3 ξ1))(*Wartość oczekiwana*)

Masa kuli

m[ξ_]:=NIntegrate[ρ[x] 4 π (x^2)/.(Subscript[ρ, c])->1,{x,near0,ξ}];
Plot[m[ξ],{ξ,0,ξ0},PlotRange->All,AxesLabel->{"ξ","m (α^3) / (Subscript[ρ, c])"},GridLines->Automatic]

Ciśnienie

P[ξ_]:=k  θ[ξ]/. sol0[[1]]/.k->1/.(Subscript[ρ, c])->1;
Plot[P[ξ],{ξ,0,ξ0},PlotRange->All,AxesLabel->{"ξ","P/(Subscript[P, c])"},GridLines->Automatic]

Dla n = 1,5

n = 1,5 w przybliżeniu opisuje strukturę Jowisza (planety gazowej), białego karła o niskiej masie, brązowego karła, a także w pełni konwekcyjny rdzeń czerwonego olbrzyma.

Clear[n,sol0,ξ0,ξ1];n=1.5;

sol0=eqLE;
sol0=NDSolve[{sol0,boundaryCond,WhenEvent[Re[θ[ξ]]==0,{ξ0=ξ,"StopIntegration"}]},θ,{ξ,start,stop}];ξ1=Re[θ'[ξ0]/. sol0[[1]]];
Plot[θ[ξ]/. sol0[[1]],{ξ,start,ξ0},PlotRange->All,AxesLabel->{"ξ","θ"},GridLines->Automatic]

Gęstość

ρ[ξ_]:=(Subscript[ρ, c]) ((θ[ξ])^n)/.sol0[[1]];
Plot[ρ[ξ]/.(Subscript[ρ, c]) ->1,{ξ,0,ξ0},PlotRange->All,AxesLabel->{"ξ","ρ /(Subscript[ρ, c])"},GridLines->Automatic]

Gęstość średnia:

⟨ρ⟩=((NIntegrate[ρ[ξ] 4 π ξ^2/.(Subscript[ρ, c])->1,{ξ,near0,ξ0}])/((4/3)π ((ξ0)^3)))(Subscript[ρ, c])

Kontrast gęstości:

((Subscript[ρ, c])/(⟨ρ⟩))

-((ξ0)/(3 ξ1))(*Wartość oczekiwana*)

Masa kuli

m[ξ_]:=NIntegrate[ρ[x] 4 π (x^2)/.(Subscript[ρ, c])->1,{x,near0,ξ}];
Plot[m[ξ],{ξ,0,ξ0},PlotRange->All,AxesLabel->{"ξ","m (α^3) / (Subscript[ρ, c])"},GridLines->Automatic]

Ciśnienie

P[ξ_]:=k  ((ρ[ξ])^(1+1/n))/.k->1/.(Subscript[ρ, c])->1;
Plot[P[ξ],{ξ,0,ξ0},PlotRange->All,AxesLabel->{"ξ","P/(Subscript[P, c])"},GridLines->Automatic]

Dla n = 3

n = 3 opisuje strukturę Słońca i gwiazd podobnych, a także w przybliżeniu rdzenie białych karłów o większej masie.

Clear[n,sol0,ξ0,ξ1];n=3;

sol0=eqLE;
sol0=NDSolve[{sol0,boundaryCond,WhenEvent[Re[θ[ξ]]==0,{ξ0=ξ,"StopIntegration"}]},θ,{ξ,start,stop}];ξ1=Re[θ'[ξ0]/. sol0[[1]]];
Plot[θ[ξ]/. sol0[[1]],{ξ,start,ξ0},PlotRange->All,AxesLabel->{"ξ","θ"},GridLines->Automatic]

Gęstość

ρ[ξ_]:=(Subscript[ρ, c]) ((θ[ξ])^n)/.sol0[[1]];
Plot[ρ[ξ]/.(Subscript[ρ, c]) ->1,{ξ,0,ξ0},PlotRange->All,AxesLabel->{"ξ","ρ /(Subscript[ρ, c])"},GridLines->Automatic]

Gęstość średnia:

⟨ρ⟩=((NIntegrate[ρ[ξ] 4 π ξ^2/.(Subscript[ρ, c])->1,{ξ,near0,ξ0}])/((4/3)π ((ξ0)^3)))(Subscript[ρ, c])

Kontrast gęstości:

((Subscript[ρ, c])/(⟨ρ⟩))

-((ξ0)/(3 ξ1))(*Wartość oczekiwana*)

Masa kuli

m[ξ_]:=NIntegrate[ρ[x] 4 π (x^2)/.(Subscript[ρ, c])->1,{x,near0,ξ}];
Plot[m[ξ],{ξ,0,ξ0},PlotRange->All,AxesLabel->{"ξ","m (α^3) / (Subscript[ρ, c])"},GridLines->Automatic]

Ciśnienie

P[ξ_]:=k  ((ρ[ξ])^(1+1/n))/.k->1/.(Subscript[ρ, c])->1;
Plot[P[ξ],{ξ,0,ξ0},PlotRange->All,AxesLabel->{"ξ","P/(Subscript[P, c])"},GridLines->Automatic]

Dla n = 5

Dla n = 5  otrzymuje się nieskończony promień, rozwiązanie odpowiada gromadom kulistym gwiazd.

Clear[n,sol0,ξ0,ξ1];n=5;

sol0=eqLE;
sol0=NDSolve[{sol0,boundaryCond,WhenEvent[Re[θ[ξ]]==0,{ξ0=ξ,"StopIntegration"}]},θ,{ξ,start,stop}];ξ0=stop;ξ1=Re[θ'[ξ0]/. sol0[[1]]];
Plot[θ[ξ]/. sol0[[1]],{ξ,start,20},PlotRange->All,AxesLabel->{"ξ","θ"},GridLines->Automatic]

Gęstość

ρ[ξ_]:=(Subscript[ρ, c]) ((θ[ξ])^n)/.sol0[[1]];
Plot[ρ[ξ]/.(Subscript[ρ, c]) ->1,{ξ,0,20},PlotRange->All,AxesLabel->{"ξ","ρ /(Subscript[ρ, c])"},GridLines->Automatic]

Gęstość średnia:

⟨ρ⟩=((NIntegrate[ρ[ξ] 4 π ξ^2/.(Subscript[ρ, c])->1,{ξ,near0,20}])/((4/3)π (20^3)))(Subscript[ρ, c]) (*Dąży do zera przez nieskończony promień*)

Kontrast gęstości:

((Subscript[ρ, c])/(⟨ρ⟩))(*Dąży do nieskończoności przez nieskończony promień*)

-((ξ0)/(3 ξ1))(*Wartość oczekiwana*)

Masa kuli

m[ξ_]:=NIntegrate[ρ[x] 4 π (x^2)/.(Subscript[ρ, c])->1,{x,near0,ξ}];
Plot[m[ξ],{ξ,0,20},PlotRange->All,AxesLabel->{"ξ","m (α^3) / (Subscript[ρ, c])"},GridLines->Automatic]

Ciśnienie

P[ξ_]:=k  ((ρ[ξ])^(1+1/n))/.k->1/.(Subscript[ρ, c])->1;
Plot[P[ξ],{ξ,0,20},PlotRange->All,AxesLabel->{"ξ","P/(Subscript[P, c])"},GridLines->Automatic]