二分法的案例
二分法+打靶法解微分方程
摘要: 先確定二分法的上下兩個端點值,這個在figure的subplot(121)中確定。
然后進行循環求解。
方程是：
diff(s,2)+2*diff(s,1)=m*s*exp(-n*r)+h*s,
where diff is a difference for r.
邊界條件:
r=0,ds/dr=0
r=R,s=5
程序：
% 二分法+打靶法解微分方程
% 方程:
% diff(s,2)+2*diff(s,1)=m*s*exp(-n*r)+h*s,
% where diff is a difference for r.
% 邊界條件:
%
% r=0,ds/dr=0
% r=R,s=5
clc;clear;close all;
m=1;
n=1;
% \copyright: zjliu
% Author's email: zjliu2001@163.com
h=1;
R=5;
fun=inline('[s(2);m*s(1)*exp(-n*r)+h*s(1)-2*s(2)/r]',...
展開 Spring-ICE 結冰算法述評-（3）水滴收集量計算
思路一：二分法
具備基礎高等數學的都知道這個方法。二分法是最基礎的數值迭代方法，通過不斷試錯二分迭代，逼近結果。
應用到我們這個問題上，就是預設一個大的水滴撞擊范圍，發射兩條水滴，然后根據撞擊結果，不斷調整發射區間直至我們設定的誤差范圍。下面這個圖就是基于該思路求解整個壁面撞擊上限水滴逼近結果。
二分法的優點就是精度高，缺點就是效率低，每個單元逐個求解過來，總共可能要計算數百條軌跡。
思路二：粗略估算法
還有個思路更簡單的辦法，就是一次計算一定水量的水滴數（比如200條），然后看每個單元打中幾條，通過比例粗略計算水滴收集量。
這個方法程序設計比較簡單，但是效率低，精度差。在實際計算時候，有大量區域沒有撞擊到水滴，水滴撞擊個數在翼面的分布呈現極大的不連續性，因此即便是按比例插值，事先還是要做濾波處理，幾次插值下來，精度根本就保證不了了。
思路三：高效高精度插值法
說到這，大家也都看出來了，水滴收集量的計算想做的漂亮，就是要同時保證計算精度和效率。那么搞出一個高精度高效率的插值方法，將變得非常關鍵。
上海交大搞了一種所謂兩級插值法，說白了，一級插值是在遠場計算少量水滴，離壁面近了以后，認為還有一個水滴和當前水滴距離比較近，通過徑向基函數插值出其遠場軌跡，在正常計算其撞擊壁面的過程。第二級插值，還是我上面說的，根據撞擊個數插值收集量，只不過還是用所謂徑向基函數。這個方法呢，怎么說呢，比較繞，寫論文會比較好看。
一般來說，我個人的研究品味還是傾向于有樸素設計思想的方法，簡單直接，又意想不到。
展開 Fortran基礎編程（入門簡介篇）
我......現如今使用的人確實很少，若不是自己對二次開發感興趣，我也不會去觸碰它~
本文先向大家簡單介紹一個二分法的fortran程序，然后推薦大家使用一些輕便的編譯器，最后推薦一下學習的書籍、論壇。
Fortran程序分為自由格式與固定格式，固定格式太太太難看了，閱讀起來也很費力，但大多數遺留下來的有限元程序是按照固定格式來編譯的（木得辦法），自由格式閱讀起來就比較友好了，接下來就用自由格式來編寫一個二分法程序，固定格式相關的將會在以后的推文中涉及。程序在文末！
本程序采用Module將代碼塊進行封裝，程序中用到的變量、function、subroutine，整理到Module中，用的時候只需要Use一下就可，如行48。程序中一些語法、注意的細節將會在以后專門設立一個Fortran的欄目中進行講解。
木木最開始用的編譯器是VS+IVF，后來發現VS占用系統資源很多，卸載也很麻煩，用起來不順手，現在我用的是Simply Fortran編譯器，小巧輕便（100M左右），很很很好用，Simply Fortran安裝包及簡單的說明手冊在文末獲取。界面大致長這個樣子：
B站有的教學視頻也用的是這個編譯器，大家有興趣可以看一下。還有一個有意思的在線編譯FORTRAN網站：http://fortran.jsrun.net/ 。界面大致長介個樣子：
接下來就是推薦書籍和論壇了。書籍毫無疑問當然是臺灣大學彭國倫先生的《Fortran 95 程序設計》，書里面有自由格式的程序，也有對應固定格式的程序，雙管齊下，大家有興趣可以點擊了解一下。
展開 《數值方法和MATLAB實現與應用》
6.3.1 二分法的分析
6.3.2 收斂準則
6.3.3 二分法的一般實現
6.4 牛頓法
6.4.1 牛頓法的收斂性
6.4.2 牛頓法的一般實現
6.5 割線法
6.6 混合法
6.7 多項式的根
6.8 小結
習題
第7章 線性代數回顧
7.1 向量
7.1.1 向量操作
7.1.2 向量的范數
7.1.3 正交向量
7.2 矩陣
7.2.1 矩陣中的每行和每列都是向量
7.2.2 對矩陣進行的操作
7.2.3 矩陣運算和向量運算的操作次數
7.2.4 矩陣的范數
7.3 向量和矩陣的數學性質
7.3.1 線性無關性
7.3.2 向量空間
7.3.3 與矩陣相關的子空間
7.3.4 矩陣的秩
7.3.5 矩陣的行列式
7.4 特殊矩陣
7.4.1 對角矩陣
7.4.2 單位矩陣
7.4.3 矩陣的逆
7.4.4 對稱矩陣
7.4.5 三對角矩陣
7.4.6 正定矩陣
7.4.7 正交矩陣
7.4.8 置換矩陣
7.5 小結
習題
第8章 解方程組
第9章 數據的最小二乘曲線擬合
第10章 插值
第11章 數值積分
第12章 常微分方程的數值積分
附錄A 特征值和特征方程組
附錄B 稀疏矩陣
參考文獻
MATLAB工具箱函數
NMM工具箱m文件函數列表
索引
展開 
提高ANSYS非線性求解收斂性能的一般方法總結
值得說明的是，當采用自動時間步長后，二分法會被自動激活。如果在一個太大時間步內收斂失敗，該特性能減半時間步，使得計算繼續。
二、Newton-Raphson 選項調整與自適應下降
在非線性計算收斂困難時，可以調整程序默認的牛頓—拉普森選項。雖然一般情況下采用采用自動默認的選項會獲得最佳的收斂特性，但也不排除會遇到其他選擇會更有效的情況。合理使用自適應下降因子也能增加某些非線性單元的收斂情況。
三、線性搜索
線性搜索可以看為是對自適應下降的一個替代，兩者不應同時使用。線性搜索一般情況下能使得分析得到收斂，但隨之帶來的是求解速度的大幅度降低，特別是針對有材料非線性的情況，一般而言，當結構采用力加載或者剛度增長的薄膜等類似分析時，可以打開線性搜索，線性搜索打開命令為LNSRCH。
四、調整收斂準則
ANSYS非線性收斂準則主要有四種，分別為力、位移、彎矩和轉角。在常用的分析中，以力作為控制加載時，可以使用殘余力的范數控制收斂；而在位移控制加載時，一般采用位移的范數控制收斂。比較常見的ANSYS非線性計算畫面如下：
FCRIT 代表的是力收斂準則；F L2是當前迭代力的收斂情況；U是位移，CRIT是收斂準則，L2是當前迭代位移收斂情況?？梢酝ㄟ^CNVTOL命令來放寬收斂準則，程序默認收斂值為0.1%,一般情況下可根據實際工程修改，例如在對鋼筋混凝土的非線性分析中，一般可放寬至5%.
五、增加迭代次數、加大荷載子步
ANSYS默認在每一荷載步計算中，最大的平衡迭代次數為25，如果在這平衡迭代次數之內不滿足收斂準則，而且自動步長是打開的，程序將使用二分法繼續計算，如果不是，將會終止分析過程，因此在實際計算中，如果發現收斂較為困難，可以適當增加平衡迭代次數的限值，以及加大荷載子步的數目。
展開 飛機結冰的那些事（2） Spring-Ice結冰軟件介紹
采用拉格朗日法，并基于二分法預測水滴收集系數。目前針對二分法計算效率低的問題，研發團隊已經提出了特征線插值法，準備用于下一版本的軟件升級。
3） 結冰模型?；诮浀銶essinger結冰模型預測結冰量。
程序界面如下圖所示：
圖1 Spring-Ice界面
實際的程序除上述三大塊原理部分還有必要的數據前后處理：
1） 前處理模塊主要針對輸入節點數目過多（高于150個點）數據進行簡化處理以提高計算效率，同時對輸入數據中的重復數據進行刪除。
2） 后處理模塊主要對產生冰形進行光順處理。
3） 數據輸出，根據輸入工況名稱輸出冰形數據。
3. 特點
和現有結冰程序相比，Spring-Ice結冰程序的特點包括以下幾個方面：
1） 計算模塊和方法盡可能與Lewice保持一致。
2） 經過冰風洞試驗的校準，該程序已內置合適的時間步，因此不需要設置時間步長等參數，只需輸入工況參數即可給出合適的結果。
3） 單個冰形的計算時間在五分鐘以內，計算效率較高。
4） 程序具有對大弦長翼型（超過1m）的準確結冰模擬能力，優于Lewice。
5） 程序具備下冰角的模擬能力。
4. 計算效果
對NACA0012、GLC-305、LTHS、NLF0414等翼型典型冰形進行計算，并與冰風洞試驗結果進行對比。相應計算結果如圖1-6所示。
展開 非線性的定義
幾何非線性：
ü 即物體在大變形時應力應變關系是非線性的，剛度連續變化，這與小變形情況下剛度是一常量不同；
ü 高速旋轉物體，離心力導致幾何形狀改變，進而產生初始應力，這種初始應力也會導致剛度的改變，通常稱為旋轉軟化；
ü 旋轉軟化效應的設置：solution>define loads >apply>structural >other>angular velocity>選中decrease stiffn選項，對應的命令為omega；
材料非線性：
ü 材料的塑性變形階段；
ü 雙線性材料
ü 多線性材料
ü 超彈性材料：多用于模擬橡膠，需采用hyper類單元；或除link，beam單元以外的其它類型的編號在180以上的單元；
ü 溫度相關的非線性材料：
ü 蠕變材料：
ü 粘塑性：
2 單元生死在狀態非線性問題中的應用
2 非線性與非彈性
非線性材料并不一定是非彈性的，即不一定是塑性的，非彈性指的是變形在撤去載荷后不可以恢復的性質；
2 大變形和大應變效應選項，nlgeom，并不是所有的非線性分析都將產生大變形；
2 非線性分析的特殊性
ü 逐步遞增載荷：非線性問題需要一系列帶校正的線性近似來求解；
ü ANSYS采用牛頓－拉普森（NR）平衡迭代方法來增加收斂性；
ü ANSYS提供了下列方法來增加問題的收斂性：自適應下降、線性搜索、自動載荷步、二分法
展開 基于ANSYS APDL去判斷一個數是不是質數 ￥15
本質ANSYS APDL迭代的過程中沒有采取什么二分法啥的，就是從2開始，到最大數-1的遍歷。
下圖為數組中的數字，因為1不為質數，數組里直接排除。
下面為APDL 經典語句運行后質數形成的數組
只是目前數組行數和前面數組一行多，無質數的由0填充。
共計273個質數。
下圖是由*vwrite輸出的質數數：
1.txt中就是運行的程序。
基于ANSYS的鋼筋混泥土復合墻板力學性能分析
圖3 網格劃分示意圖 （a）solid65單元（b）beam188單元
三、載荷的加載
對于ANSYS求混泥土的極限載荷有兩種加載方式：
力加載：通過對應的方法（比如特征值屈曲）估計結構的極限載荷的大致范圍，然后給結構施加一個稍大的載荷，打開自動載荷步二分法進行非線性靜力分析，最后計算會因不收斂而終止，則倒數第二個子步對應的就是結構的極限載荷；另外也可以采用弧長法，采用足夠的步數（弧長法一直可以分析到極限載荷承載力之后的過程）同樣可以從繪制的載荷位移曲線或計算結果中找出結構的極限載荷。
位移加載：給結構施加一個大的位移，打開自動載荷步二分法進行非線性分析，保證足夠的子步數，這樣也可以通過繪制載荷位移曲線或者查看相應結果文件知道結構的極限載荷。
本文采用位移加載的方式。相對于力加載，位移加載在求解速度上更有優勢。原因如下：涉及到塑性求解時，有材料的應力——應變曲線可知，材料的彈性模量不是定值，每一步求解中的材料彈性模量利用的是上一步的彈性模量值，而加載力之后，根據上一次的彈性模量求位移，進而求得應變，根據材料的應力-應變曲線，得到材料此位移下材料的真實彈性模型，與上一步的彈性模量指進行對比，進而更新第一次用的彈性模量，反復更新才能計算收斂。而位移加載，直接算得應變，根據材料的應力-應變曲線，可得到彈性模量，進而得到力，不用反復更新彈性模量，求解效率大大提高。
本文采用第二種求解方式，在復合墻板的正中加一個大的位移，根據底部固定節點的反力從而得到加載復合墻板上的載荷，同時利用acel命令流對鋼筋混泥土的重量進行加載。相應的加載曲線如下圖4所示。
圖4 墻板的位移加載示意圖
四、求解器的控制
ANSYS混泥土的非線性計算收斂是非常困難的。主要影響因素有網格密度、子步數、收斂準則等。
網格密度：網格密度適當能夠收斂。
展開 基于ANSYS的鋼筋混泥土復合墻板力學性能分析
圖3 網格劃分示意圖 （a）solid65單元（b）beam188單元
三、載荷的加載
對于ANSYS求混泥土的極限載荷有兩種加載方式：
力加載：通過對應的方法（比如特征值屈曲）估計結構的極限載荷的大致范圍，然后給結構施加一個稍大的載荷，打開自動載荷步二分法進行非線性靜力分析，最后計算會因不收斂而終止，則倒數第二個子步對應的就是結構的極限載荷；另外也可以采用弧長法，采用足夠的步數（弧長法一直可以分析到極限載荷承載力之后的過程）同樣可以從繪制的載荷位移曲線或計算結果中找出結構的極限載荷。
位移加載：給結構施加一個大的位移，打開自動載荷步二分法進行非線性分析，保證足夠的子步數，這樣也可以通過繪制載荷位移曲線或者查看相應結果文件知道結構的極限載荷。
本文采用位移加載的方式。相對于力加載，位移加載在求解速度上更有優勢。
原因如下：涉及到塑性求解時，有材料的應力——應變曲線可知，材料的彈性模量不是定值，每一步求解中的材料彈性模量利用的是上一步的彈性模量值，而加載力之后，根據上一次的彈性模量求位移，進而求得應變，根據材料的應力-應變曲線，得到材料此位移下材料的真實彈性模型，與上一步的彈性模量指進行對比，進而更新第一次用的彈性模量，反復更新才能計算收斂。
而位移加載，直接算得應變，根據材料的應力-應變曲線，可得到彈性模量，進而得到力，不用反復更新彈性模量，求解效率大大提高。
本文采用第二種求解方式，在復合墻板的正中加一個大的位移，根據底部固定節點的反力從而得到加載復合墻板上的載荷，同時利用acel命令流對鋼筋混泥土的重量進行加載。相應的加載曲線如下圖4所示。
圖4 墻板的位移加載示意圖
四、求解器的控制
ANSYS混泥土的非線性計算收斂是非常困難的。主要影響因素有網格密度、子步數、收斂準則等。
網格密度：網格密度適當能夠收斂。
展開 基于基于注塑模CAD/CAE的自由曲面的裁剪的自由曲面的裁剪
對裁剪曲面中每個Patch在求交完畢后隨即將其分割成2個或多個Patch，不論是2個或多個Patch，必須按一定規則將其分為2類，一類保留，另一類裁剪掉，因而稱為及時分片法.如圖5所示，一般情況下一個Patch被一分為二，則這2個新生成的Patch自然分為上述2類.但在某些情況下，也可能被分為多個Patch，如圖6所示.用逐次二分法將其分解為多個Patch后，本文采用交線包容法，同樣可將這多個Patch分為2類，如圖7所示.所謂交線包容法是將包含原Patch內所有交線的多邊形作為一類，如圖7(a)所示，其它的作為另一類，如圖7(b)所示.
圖 5 1個平面片裁剪后生成2個平面片
圖 6 逐次二分法分解平面片
圖 7 交線包容法分類
步驟2. 將新生成的Patch加入到Patch鏈表中，去掉原Patch.
步驟3. 處理如圖8所示的特殊情況.一個Patch中出現內環，無法用2個多邊形來描述分割后的圖形，對這種情況本文采取加密離散網格的方法，避免出現內環.
圖 8 平面片內環情況
5、生成裁剪曲面的自動搜索算法
圖9為求交后的結果，被裁剪曲面中有一部分Patch是將原Patch裁剪生成的，稱為新Patch，如Patch1，2，3，而有一部分是沒被裁剪到的，稱為原Patch，如Patch4，5，下面分別討論這2種Patch的搜索算法.
圖 9 曲面求交后的搜索
5.1 新Patch的搜索算法
步驟1. 由一個新Patch找到與之具有相鄰公共邊的新Patch，如圖9中由Patch3可找到Patch1，2.
步驟2.
展開 
基于 MIST 方法提出了點陣結構的尺寸優化算法（用于制造增材優化等課題） ￥800
水平平面對應的目標值 t 取決于體積約束，如果當前迭代步的實體材料過多則增加 t 的數值使水平面向上移動，反之則降低 t 的數值使水平面向下移動，t 在每一步迭代步中的具體數值可以通過二分法等方式計算得到。與 SIMP 方法中的密度類似，MIST 方法定義了一種體積權重值來描述材料種類，體積權重為 1 表示實體材料，提及權重為 0 表示孔洞材料。體積權重更新準則如下，當單元對應的節點的響應值全都高于水平平面對應的響應值時讓單元朝著實體方向演變（體積權重變向 1），當單元對應節點的響應值全都小于水平平面對應的目標值時讓單元朝著孔洞方向演變（體積權重變向 0），當單元的部分節點的響應值小于水平平面對應的目標值而另一部分大于時，則通過響應面投影面積的方法計算單元的貢獻量。
優化問題的目標是整體結構的剛度最大，實現方法是使整體結構的總應變能最小，約束包括體積約束和設計變量的上下限約束。結構桿件的初始橫截面是通過體積約束確定的。算法的優化列式如下：
圓環優化的命令流如下（收費內容，也可以聯系qq：290539594）.
展開 ANSYS Workbench掃地機器人邊刷靜力分析【原創】
網格劃分
首先劃分較粗糙網格進行試算，檢驗材料參數、零件裝配、接觸對以及載荷等是否正常，再采用二分法調整單元尺寸或對局部進行細化，直至結果滿足精度要求。
求解設置
接觸問題不收斂可能有眾多原因，如材料、接觸對的定義、網格尺寸等因素有關，但更多時候是求解設置不當導致分析不收斂。建議打開大變形開關，采用自動時間步，適當調整初始/最小/最大求解子步以及接觸剛度，通常情況下即可取得收斂結果。
分析結果
刷頭與地面的接觸動畫如下圖
應力分布
接觸力的變化
刷頭的位移響應
最開始沒有找到Workbench中的解決方法，只好將結果又導入APDL中，將刷頭的位移按求解子步寫出，再繪成曲線。
后面發現其實Workbench中處理更方便，直接定義一個User Defined Result，省去編寫APDL代碼寫出結果的麻煩。
3.總結
刷頭與地面的靜力接觸力分析，是一個基礎、典型的接觸分析，涉及局部坐標系、遠程點/遠程位移、接觸對定義、Named Selections、非線性求解設置，以及后處理中變形、應力、接觸力、節點位移響應等應用。
展開 今日Nature：電子摻雜氧化銅超導體中的三維空間電荷激發
這種二分法強調了面外電荷動力學的重要性。目前研究人員已經發現，在光學可達到的有限動量范圍內，面外電荷動力學在正常狀態下是非相干的。
【成果簡介】
今日，來自斯坦福大學的Z. X. Shen, T. P. Devereaux,和W. S. Lee（共同通訊）聯合在Nature發表文章，題為“Three-dimensional collective charge excitations in electron-doped copper oxide superconductors”。作者使用共振非彈性X射線散射來探索布里淵區所有三維空間的電荷動態。電子摻雜氧化銅激發(模式)的極化分析揭示了它們的電荷來源。激發沿著面內和面外兩個方向擴散，這揭示了它的三維性質。面外分散的周期性對應于相鄰CuO2平面之間的間距，而不是晶體c軸晶格常數，這表明面外庫侖相互作用是相干面外電荷動力學的原因。作者觀察到的特性是渴望已久的“聲等離子體激元”的特征，這是分層系統預測的不同電荷集合模式的一個分支，并被認為在調解高溫超導性方面將發揮重要作用。
【圖文簡介】
圖1. 層狀電子氣中的等離子體激元以及電子摻雜氧化銅的電荷激發
圖2. 區域中心的激發三維性
圖3.
展開 ANSYS中弧長法的原理
為了搞清楚得到的最大荷載是否是結構真實的極限荷載，需要用弧長法來幫幫忙：
用弧長法進行預分析，得到結構屈曲荷載近似值(預測數值)，再用傳統的二分法計算，兩個結果進行比對看是否一致；
使用弧長法計算，計算中手動修改弧長半徑，再看結果的變異情況。
繪制出結構的荷載-位移曲線，探討曲線變化點的原因，從而確定數據是否可信。
于是，弧長法非常值得研究。
何為弧長法？
弧長法是一種非線性求解的迭代控制方法，由于其可以解決在荷載和位移增量均不確定的情況下，生成變化的增量值，并能很好地追蹤結構加載路徑而具有很高的“聲望”。關于弧長法的原理，推薦參考《非線性分析弧長法的讀書報告》、陸新征老師學生時代的作業：《基于預處理技術和弧長法的非線性方程通用求解子程序總結報告 》，以及Yusd的博文《弧長法（Riks Method）的基本原理》。喜歡編程的話，還可以參考他的另一篇文章《弧長法(Riks method)通用求解程序》。英文資料可閱讀蘇黎世聯邦理工學院結構工程研究所Prof. Dr. Eleni Chatzi的講稿：《The Finite Element Method for the Analysis of Non-Linear and Dynamic Systems》
在這里只強調一下弧長法的一個獨特的優勢：
見下圖，當微小荷載增量可以引起顯著的位移突躍，則成為荷載控制的急速通過(Snap though)；當微小位移增量可以產生顯著的荷載突躍，則成為位移控制的急速返回(snap back)——這種復雜的加載路徑，單純荷載或位移控制已經難以奏效，而弧長法則可輕松處理此類問題。
轉自公眾號——ANSYS學習與應用
旨在分享，若侵即刪.
展開 