實驗目的
用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 %二分法改良 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
運行示例(同樣沒有控制輸出)
明顯地,迭代次數減小許多。
繼續優化代碼,以獲得區間里所有的根(關鍵是對有根區間的判斷與記憶)

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
運行示例(仍然沒有控制輸出,筆者忽然覺得自己中了月亮的毒!!!)
使用二分法解得在[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!!!待解決——有空再說。)