Matlab非線性方程數值解法


實驗目的

Matlab實現非線性方程的二分法、不動點迭代法

實驗要求

1. 給出二分法算法和不動點迭代算法

2. Matlab實現二分法

3. 用Matlab實現不動點迭代法

實驗內容

(1)在區間[0,1]上用二分法和不動點迭代法求的根到小數點后六位。

2)二分法的基本思想:逐步二分區間[a,b],通過判斷兩端點函數值的符號,進一步縮小有限區間,將有根區間的長度縮小到充分小,從而,求得滿足精度要求的根的近似值。

3)不動點迭代法基本思想:已知一個近似根,構造一個遞推關系(迭代格式),使用這個迭代格式反復校正根的近似值,計算出方程的一個根的近似值序列,使之逐步精確法,直到滿足精度要求(該序列收斂於方程的根)。

實驗步驟

(1)二分法算法與MATLAB程序(二分法的依據是根的存在性定理,更深地說是介值定理)。

 

 

 

MATLAB程序,

 1 %二分法
 2 %輸入:f(x)=0的f(x),[a,b]的a,b,精度ep
 3 %輸出:近似根root,迭代次數k
 4 function [root,k]=bisect(fun,a,b,ep)
 5 if nargin>3
 6     elseif nargin<4
 7         ep=1e-5;%默認精度
 8     else
 9         error('輸入參數不足');%輸入參數必須包括f(x)和[a,b]
10 end
11 if fun(a)*fun(b)>0%輸入的區間要求
12     root=[fun(a),fun(b)];
13     k=0;
14     return;
15 end
16 k=1;
17 while abs(b-a)/2>ep%精度要求
18     mid=(a+b)/2;%中點
19     if fun(a)*fun(mid)<0
20         b=mid;
21     elseif fun(a)*fun(mid)>0
22         a=mid;
23     else
24         a=mid;b=mid;
25     end
26     k=k+1;
27 end
28 root=(a+b)/2;
29 end
二分法1

運行示例(並未對輸出格式做控制,由於精度要求,事后有必要控制輸出的精度):

 優化代碼,減小迭代次數(在迭代前,先搜尋更適合的有根區間)

 1 %二分法改良
 2 %在一開始給定的區間中尋找更小的有根區間
 3 %輸入:f(x)=0的f(x),[a,b]的a,b,精度ep
 4 %輸出:近似根root,迭代次數k
 5 %得到的根是優化區間里的最大根
 6 function [root,k]=bisect3(fun,a,b,ep)
 7 if nargin>3
 8     elseif nargin<4
 9         ep=1e-5;%默認精度
10     else
11         error('輸入參數不足');%輸入參數必須包括f(x)和[a,b]
12 end
13 %定義划分區間的分數
14 divQJ=1000;
15 %等分區間
16 tX=linspace(a,b,divQJ);
17 %計算函數值
18 tY=fun(tX);
19 %找到函數值的正負變化的位置
20 locM=find(tY<0);
21 locP=find(tY>0);
22 %定義新區間
23 if tY(1)<0
24 a=tX(locM(end));
25 b=tX(locP(1));
26 else
27 a=tX(locP(end));
28 b=tX(locM(1));
29 end
30 if fun(a)*fun(b)>0%輸入的區間要求
31     root=[fun(a),fun(b)];
32     k=0;
33     return;
34 end
35 k=1;
36 while abs(b-a)/2>ep%精度要求
37     mid=(a+b)/2;%中點
38     if fun(a)*fun(mid)<0
39         b=mid;
40     elseif fun(a)*fun(mid)>0
41         a=mid;
42     else
43         a=mid;b=mid;
44     end
45     k=k+1;
46 end
47 root=(a+b)/2;
48 end
二分法2

運行示例(同樣沒有控制輸出)

明顯地,迭代次數減小許多。

繼續優化代碼,以獲得區間里所有的根(關鍵是對有根區間的判斷與記憶)

  1 %二分法改良
  2 %在一開始給定的區間中尋找更小的有根區間
  3 %輸入:f(x)=0的f(x),[a,b]的a,b,精度ep
  4 %輸出:近似根root,迭代次數k
  5 %例子:
  6 %fun=inline('x.^3-x.^2-4*x+4');a=-3;b=3;ep=1e-5;
  7 %fun=inline('x.^2-1');a=-2;b=2;ep=1e-5;
  8 %fun=inline('-x.^2+1');a=-2;b=2;ep=1e-5;
  9 %fun=inline('x.^2-3*x+2');a=0;b=3;ep=1e-5;
 10 function [root,k,x,y]=bisect4(fun,a,b,ep)
 11 if nargin>3
 12     elseif nargin<4
 13         ep=1e-5;%默認精度
 14     else
 15         error('輸入參數不足');%輸入參數必須包括f(x)和[a,b]
 16 end
 17 %定義划分區間的分數
 18 divQJ=100;
 19 %等分區間
 20 tX=linspace(a,b,divQJ);
 21 %計算函數值
 22 tY=fun(tX);
 23 %找到函數值的正負變化的位置
 24 locM=find(tY<0);
 25 locP=find(tY>0);
 26 %得到的符號序列,找其中的間斷點
 27 z=0;
 28 k=1;
 29 for i=1:length(locM)-1%掃描-符號序列
 30     %滾筒算法
 31     if(i==1 && locM(1)~=1)
 32         z=z+1;
 33         x(z)=locM(i);
 34     end
 35     tmp=locM(i+1)-locM(i);
 36     if(tmp~=1)
 37         z=z+1;
 38         x(z)=locM(i);
 39         z=z+1;
 40         x(z)=locM(i+1);
 41     end
 42 end
 43 
 44 l=0;
 45 for i=1:length(locP)-1%掃描+符號序列
 46     %滾筒算法
 47     if(i==1 && locP(1)~=1)
 48         l=l+1;
 49         y(l)=locP(i);
 50     end
 51     tmp=locP(i+1)-locP(i);
 52     if(tmp~=1)
 53         l=l+1;
 54         y(l)=locP(i);
 55         l=l+1;
 56         y(l)=locP(i+1);
 57     end
 58 end
 59 if(z==l-1)
 60     z=z+1;
 61     x(z)=locM(end);
 62 elseif(z==l+1)
 63     l=l+1;
 64     y(l)=locP(i);
 65 end
 66 
 67 if(z==0 && l==0)
 68     if(locM(end)<locP(1))
 69         x(1)=locM(end);
 70         y(1)=locP(1);
 71     else
 72         x(1)=locM(1);
 73         y(1)=locP(end);
 74     end
 75     z=1;
 76 elseif(z==0)
 77     x(1)=locM(1);
 78     x(2)=locM(end);
 79     z=2;
 80 elseif(l==0)
 81     y(1)=locP(1);
 82     y(2)=locP(end);
 83     z=2;
 84 end
 85 for i=1:z
 86 a=tX(x(i));
 87 b=tX(y(i));
 88 if fun(a)*fun(b)>0%輸入的區間要求
 89     root=[fun(a),fun(b)];
 90     k=0;
 91     return;
 92 end
 93 if(a>b)
 94     tmp=b;
 95     b=a;
 96     a=tmp;
 97 end
 98 while abs(b-a)/2>ep%精度要求
 99     mid=(a+b)/2;%中點
100     if fun(a)*fun(mid)<0
101         b=mid;
102     elseif fun(a)*fun(mid)>0
103         a=mid;
104     else
105         a=mid;b=mid;
106     end
107     k=k+1;
108 end
109 root(i)=(a+b)/2;
110 end
111 end
二分法3

運行示例(仍然沒有控制輸出,筆者忽然覺得自己中了月亮的毒!!!)

 

 

 

 

 使用二分法解得[0,1]上的根為0.090526。(終於在最后回答問題時做了<_>     控制輸出的方法:在調用函數前,在控制台輸入format long;得到一long型的數據,根據精度讀取。)

(2)不動點迭代算法與MATLAB程序。

 

 1 %不動點迭代法
 2 %輸入:fun函數,初始值x0,精度ep,最大迭代次數Nmax
 3 %輸出:近似根root,迭代次數k
 4 function [root, k]=iterate(fun,x0,ep,Nmax)
 5 if nargin<4
 6     Nmax=500;
 7 end
 8 if nargin<3
 9     ep=1e-5;
10 end
11 x=x0;
12 x0=x+2*ep;
13 k=0;
14 while abs(x0-x)>ep && k<Nmax
15     x0=x;
16     x=fun(x0);
17     k=k+1;
18 end
19 root=x;
20 if k==Nmax
21     warning('已達到迭代次數上限');
22 end
23 end
不動點迭代法

運行示例(對迭代方程的簡單推導得到:):

 就沒優化的二分法與不動點迭代法比較,明顯地迭代法的迭代次數更少。

使用不動點迭代法,求得近似根0.090525

 小結

就上述兩個程序而言,二分法的編程更具挑戰。在優化二分法的程序過程中,發現二分法的關鍵就在於對所給定區間的討論划分。得到一個能求解給定區間里的所有零點的程序的關鍵是(分類討論的思想),對符號變化的記錄與使用,換言之,解決這個問題的關鍵就在於分支結構的使用。

(不可否認,第三個二分法程序仍然有bug!!!待解決——有空再說。)


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM