用數值模擬行星運動,很難得到穩定解?

時間 2021-05-29 22:57:51

1樓:qfzklm

我在這裡貼乙個我之前遇到的相關聯的問題。。

為什麼這個演算法誤差的看起來這麼小? - 知乎

簡而言之,確實需要專門做理論推導才能設計出來合適的差分方程來進行計算機模擬。

2樓:Joseph

不能忍受忽略了所有物理常數的物理公式!!

應該是這樣:

import

numpy

asnp

import

matplotlib.pyplot

asplt

defrun(N

):t=np

.arange

(0.0,N

,0.1)x

=np.zeros(N

*10)=

np.zeros(N

*10)ax

=np.zeros(N

*10)y

=np.zeros(N

*10)vy

=np.zeros(N

*10)ay

=np.zeros(N

*10)r

=np.zeros(N

*10)G

=1m=

1d=1

defAcceleration(d

,G,m

):return-G

*m/(

d**2)

x[0]

=0.5*d

[0]=

-0.2*d

y[0]

=0.0*d

vy[0]

=1.63*d

r[0]

=np.sqrt(x

[0]**

2+y[

0]**2

)ax[0

]=Acceleration(r

[0],G

,m)*

x[0]

/r[0

]ay[0

]=Acceleration(r

[0],G

,m)*

y[0]

/r[0

]delta

=0.01/(

np.sqrt(ax

[0]**

2+ay[

0]**2

))foriin

range(1

,N*10

):[i]

=[i-

1]+ax

[i-1

]*deltavy[

i]=vy

[i-1

]+ay[

i-1]

*deltax[

i]=x

[i-1

]+[i

-1]*

deltay[

i]=y

[i-1

]+vy[

i-1]

*deltar[

i]=np

.sqrt(x

[i]**

2+y[

i]**2

)ax[i

]=Acceleration(r

[i],G

,m)*

x[i]

/r[i

]ay[i

]=Acceleration(r

[i],G

,m)*

y[i]

/r[i

]plt

.title

("Simulation Times: %d"

%int(N

*10))plt

.plot(x

,y,'o'

)plt

.show

()if

__name__

=="__main__"

:run

(2500

)作個圖發現軌跡跑遠了:

說明步長設定為0.01太大了。改小點試試!

delta

=0.001/(

np.sqrt(ax

[0]**

2+ay[

0]**2

))果然好多了:收工。

3樓:Comzyh

我記得解微分方程,第一講就是如果搞不定就麼上計算機。

直接模擬有系統誤差,比如你就算模擬乙個二次函式,如果直接按定義來模擬,由於步長不是0,一定會偏向外側/內側。

估計用容格庫塔會好很多,一味減小步長不是正路。

4樓:

但是很慚愧沒有弄清穩態的本質,也就是回答偏了

題主問題的根本原因,應該是 @浪客所說的「 長期模擬行星運動要用專門設計的保能量的辛演算法常規的演算法會導致能量累積或衰減,高階龍格庫塔或者隱式方法只是誤差累積的慢一些。」

我專門查了一點相關文獻,的確如此。

題主問題中的穩定性,應該是指的「不出現隨著時間一直累積的數值計算誤差」,這個跟微分方程數值求解過程本身的穩定性,是兩個不同的概念。

這裡也非常感謝@ Sun AO,@KL Hahns等知友對我答案的指正。

下面是原答案,請大家隨便看看就好——因為答偏了。

原答案

個人也認為主要是問題出在所選用的數值演算法上,比如,時間步長不夠小。

隱式演算法對數值求解穩定性方面很有優勢,也可以控制求解精度,同時,可以進行變步長計算(在要求的精度下),以節約求解時間

如上面的知友所述,可以考慮隱式尤拉法。隱式尤拉法是最簡單的隱式常微分方程(組)演算法,但需要利用牛頓法等迭代求解非線性方程(組)。題主有興趣可以搜搜相關文獻和程式。

數值計算方法的重要性,其實遠超大部分科學和工程技術人員的想象

原答案

5樓:

兩體問題一般是二次曲線,在遠離拋物線、圓、直線這些臨界情形時不涉及初始狀態。是橢圓就是橢圓,是雙曲線就是雙曲線。

你用的是顯式尤拉折線法。

隱式尤拉折線法無條件穩定。

顯式尤拉折線法條件穩定,且需要步長充分小。

不過這個穩定性指的是小的初始誤差帶來的誤差的,而你口中的穩定性似乎只是捨入誤差和截斷誤差。

不過,這一步的捨入誤差和截斷誤差放到下一步就是初始的誤差,所以套用穩定性的分析也沒錯。

改成極座標,並用龍格庫塔四階方法去計算,會好些。

如何選擇數值模擬軟體?

楊胖胖 用FLUENT做過類似煤自燃問題,組分 反應產熱 蒸發吸熱源項都是用UDF定義。水蒸發較為麻煩,一般認為和分壓以及溫度耦合,我在計算時直接設定乙個蒸發速率計算的。如果有高明做法,還請指教。 z1840393 你這個問題用不著軟體。多孔介質輸運,涉及多個過程,每個過程估計你都有自己的方程,而這...

每次模擬行測都40多分 怎麼學都是這樣

樹霧 你這個情況還是基礎知識點掌握的不夠紮實啊!不然分數應該在五六十分的,然後在此基礎上學習做題技巧和做題方法,基礎知識 多練習 多總結,行測突破瓶頸不是問題。像我,我的言語部分當時學的很差,但是我們省這塊內容佔比又很大,第一遍跟著粉筆系統班學完基礎,開始做題的時候,錯的很多,然後把每種題型的錯題整...

作為一名用fluent做數值模擬的研究生,可為就業做從哪些準備?

丫頭 題主研二,正是學完基礎課程開始充實自己的時間,如果不想繼續深造讀博士,那我覺得乙個軟體倒是也可以了,建議如下 了解所學學科的基礎理論課程,想必題主專業是能源之類的吧,風輪機是不是與流體相關的呢,所以注重一下這些基礎的補寄。因為我是涉水專業,所以在日常的模擬中,除了fluent這種軟體,觸類旁通...