台北總公司／台北縣汐止市新台五路一段 112 號 10 樓 A 棟

1. 書名 Java於資料結構與演算法之實習應用 書號　PG20098 原著　河西朝雄著 譯者　周明憲/徐堯譯 台北總公司／台北縣汐止市新台五路一段112號10樓A棟 Building A, 10F, 112, Shin-Tai-Wu Road Sec. 1, Shijr, Taipei, Taiwan. 電話／02-26962869　傳真／02-26962867 台中辦事處／台中市文心路一段540號4樓-1 4-1F, 540, Wen Shing Road, Sec. 1, Taichung, Taiwan. 電話／04-23287870　傳真／04-23108168 博碩網址：http://www.drmaster.com.tw

2. 第二章 學習重點 • 電腦（電子計算機）原本的開發目的即為數值分析與運算。因此數值分析的演算法從一開始就有相關的研究，而其理論也有完整的體系。 • 本章並不是在解析函數的定積分，而是在說明數值積分法、泰勒展開式、非線性方程式與聯立方程式的解法、內插、最小平方法、多位數計算與作業研究(OR, Operation Research)領域的線性規劃法

3. 2-0何謂數值運算 • 「數值運算」狹義來說，是指方程式的解法、數值積分法之類的「數值分析」，廣義來說，則包括了「統計分析」、「作業研究」。

4. 非線性方程式的解法 2分法*、牛頓法* 方程式的解法 聯立1次方程式的解法 　高斯‧喬丹法*、高斯消去法* 矩陣運算 　高斯‧塞德爾法 微分方程式的解法 數值積分 尤拉法、倫格‧庫克法 數值分析 　梯形法則*、辛普森法則*、高斯公式 內插 拉格朗基內插*、牛頓內插*、史普蘭恩內插 初等函數計算 泰勒展開式* 高速傅立耶變換 基礎統計 標準偏差、相關 估計、檢測 統計分析 檢測* 遞迴分析 再遞迴分析、最小平方法* 多變量分析 主成份分析 規劃法 線性規劃法(LP) *、2次規劃法 作業研究 模擬 佇列 系統應用 PERT、最短路徑問題

5. 2-1 亂數 • 亂數(random number)是指沒有任何規則可循而隨機產生的數。對於由某一公式運算所求得的亂數則稱為假隨機亂數(pseudo-random number)。由計算所產生的亂數數列可能會有產生同一數字的時候，因而發生反覆計算情形，這種反覆計算的循環如果非常大，在實務上會將其視為亂數，我們稱之為「假隨機亂數」

6. 2-1 線性同餘法 • 先由適當的初始值開始，然後以下面的式子 xn=(Axn-1+C)%M • 依序產生0～M範圍的值。A，C，M為適當的整數，%為求餘的演算子。以除以M，8除以A餘5的數、C為奇數的條件，在週期M中依次產生0～M-1的整數。像這種各值皆平均（同一出現頻率）出現的亂數稱為一致亂數。在下列的程式中，設：A=109,C=1021,M=32768,x0=13。

7. 2-1 線性同餘法 class Rei08Panel extends Panel { int rndnum=13; // 亂數的初始值 int irnd() { // 0～32767的整數亂數 rndnum=(rndnum*109+1021) % 32768; return rndnum; } public void paint(Graphics g) { int i,y=1; g.drawString("所求得的一致亂數",0,20); for (i=0;i<100;i++) { if (i%8==0) y++; g.drawString(""+irnd(),(i%8)*60,y*20); } } }

8. 2-2 數值積分 • 例題9 梯形法則(Trapezoidal Rule)定積分 : 以梯形法 則求出函數的定積分 • 不以分析的手法（數學教科書的方法）求出函數的定積分數式，而是將其分割成微小區間，再求出近似值的方法稱為「數值積分」。

9. 2-2 數值積分 • 下列計算式為以梯形法則求 的方法。 • 圖中將a,b區間分成n個梯形，各梯形的面積合計如下：

10. 2-2 數值積分 public void actionPerformed(ActionEvent e) { int k; double a,b,n,h,x,s,sum; a=Double.parseDouble(tf1.getText()); b=Double.parseDouble(tf2.getText()); n=50; // a～b間的分割數 h=(b-a)/n; // 區間幅度 x=a;s=0; for (k=1;k<=n-1;k++) { x=x+h; s=s+f(x); } sum=h*((f(a)+f(b))/2+s); String result=" /"+b+"\n"+ " | sqrt(4-x*x) ="+sum+"\n"+ " /"+a; ta.setText(result); }

11. 2-2 辛普森法則定積分 • 以辛普森法則求出函數的定積分 • 梯形法則是用直線使 的微小區間近似起來，而辛普森法則則是用2次曲線使此區間近似起來。

12. 2-2 辛普森法則定積分 • 通過 3點的2次曲線方程式 如下： • 這個 的 間的積分值在數學上可寫成：

13. 不過， ，將此式用至a～b的區間，則 而

14. 2-2 辛普森法則定積分 • 在撰寫程式時，奇數的合計項 f0與偶數項 fe的合計要分開計算，可用下列式子計算： f(b-h)為奇數項的1個多餘項(x2n-1)。

15. 2-2 辛普森法則定積分 public void actionPerformed(ActionEvent e) { int k; double a,b,n,h,fo,fe,sum; a=Double.parseDouble(tf1.getText()); b=Double.parseDouble(tf2.getText()); n=50; // a～b之間的分割數 h=(b-a)/(2*n); // 區間幅度 fo=0;fe=0; for (k=1;k<=2*n-3;k=k+2) { fo=fo+f(a+h*k); fe=fe+f(a+h*(k+1)); } sum=(f(a)+f(b)+4*(fo+f(b-h))+2*fe)*h/3; String result=" /"+b+"\n"+ " | sqrt(4-x*x) ="+sum+"\n"+ " /"+a; ta.setText(result); }

16. d s 2-3 泰勒展開式 • 以泰勒(Taylor)展開式計算 ex可得到下列的式子： • 上列的式子由於是無限級數，因此在實際的計算上須於有限次數上截斷。截斷的條件為將項之前的和設為d，項之前的和設為s時，得到 • 的結果的時候。 稱為截斷誤差， 稱為相對截斷誤差。EPS=1e-8值則視所需的精確度來設定。若，則精確可設為8位數左右。

17. double myexp(double x) { final double EPS=1e-08; double s=1.0,e=1.0,d=1.0; int k; for (k=1;k<=200;k++) { d=s; e=e*x/k; s=s+e; if (Math.abs(s-d)<EPS*Math.abs(d)) // 截斷誤差 return s; } return 0.0; // 沒有收斂時 } public void paint(Graphics g) { double x; g.drawString("xmyexp(x)exp(x)",0,20); for (x=0;x<=40;x=x+10) { g.drawString(""+x,20,(int)x*2+40); g.drawString(""+myexp(x),80,(int)x*2+40); g.drawString(""+Math.exp(x),260,(int)x*2+40); } }

18. 2-3 負值版ex • 一般而言，在泰勒展開式內，展開時如果越靠近中心，便越能賦予好的近似值，但如離中心太遠，則誤差就會變得越大。 • 尤其是ex內的x值為負時，會對真值反覆進行大量的加算與減算，發生位數掉落的情形，使得誤差愈加擴大。 • 以例題10的程式計算e-40，以122項收斂，則 • 而得到與真值4.2485354 x10-18完全不同的值，因此x如為負的話，則以 • 求出。

19. double myexp(double x) { final double EPS=1e-08; double s=1.0,e=1.0,d=1.0,a; int k; a=Math.abs(x); for (k=1;k<=200;k++) { d=s; e=e*a/k; s=s+e; if (Math.abs(s-d)<EPS*Math.abs(d)) { // 截斷誤差 if (x>0) return s; else return 1.0/s; } } return 0.0; // 沒有收斂時 } public void paint(Graphics g) { double x; g.drawString("xmyexp(x)exp(x)",0,20); for (x=-40;x<=40;x=x+10) { g.drawString(""+x,20,(int)(x+40)*2+40); g.drawString(""+myexp(x),80,(int)(x+40)*2+40); g.drawString(""+Math.exp(x),260,(int)(x+40)*2+40); } }

20. 2-3 以泰勒展開式求cos x • 在泰勒展開式內，展開越靠近中心，便越能賦予好的近似值，因此x的值必須在0～2π的範圍內修整計算。而mycos的引數則設定為度。

21. class Dr10_2Panel extends Panel { double mycos(double x) { final double EPS=1e-08; double s=1.0,e=1.0,d=1.0; int k; x=(x-(int)x+(int)x%360)*3.1415927/180; // x的值須收斂在0～2π內 for (k=1;k<=200;k=k+2) { d=s; e=-e*x*x/(k*(k+1)); s=s+e; if (Math.abs(s-d)<EPS*Math.abs(d)) // 截斷誤差 return s; } return 9999.0; // 沒有收斂時 } public void paint(Graphics g) { double x; g.drawString("xmycos(x)cos(x)",0,20); for (x=0;x<=180;x=x+10) { g.drawString(""+x,20,(int)x*2+40); g.drawString(""+mycos(x),80,(int)x*2+40); g.drawString(""+Math.cos(x*3.1415927/180),260,(int)x*2+40); } }

22. 2-4 非線性方程式的解法 • 例題11 2分法以2分法求出方程式的根 • 1次方程式（亦即圖形上的直線）以外的方程式稱為非線性方程式。求這種方程式的根的方法中，有1個方法為2分法。

23. 2-4 非線性方程式的解法 1) 將根的左右2點a,b之初始值設為low與high。 2) 以求low與high的中點x。 3) 若，則根位於x的左邊，high=x，並將上限縮小至一半若，則根位於x的右邊，low=x，並將下限縮小至一半。 4) 設為求0或時的x值的根，否則反覆執行(2)以後各項步驟。 也就是說，2分法將資料範圍分割成2半，再反覆搜尋根位在哪半邊，然後逐步將搜尋範圍朝根的位置縮小。

24. class Rei11Panel extends Panel { double f(double x) { return x*x*x-x+1; } public void paint(Graphics g) { final double EPS=1e-08; // 截斷誤差 final int LIMIT=50; // 截斷次數 double low,high,x; int k=1; low=-2.0;high=2.0; for (k=1;k<=LIMIT;k++) { x=(low+high)/2; if (f(x)>0) high=x; else low=x; if (f(x)==0 || Math.abs(high-low)<Math.abs(low)*EPS) { g.drawString("x="+x,10,20); break; } } if (k>LIMIT) g.drawString("沒有收斂",10,20); } }

25. 2-4 以牛頓法求出方程式的根 (1)將根附近的適當的值x0設為初始值。 (2)畫出y=f(x)的x=x0的接線，將其與x軸的交點設為x1，然後以同樣的步驟求出x2, x3, ..., xn-1, xn。 (3)將 時的xn的值設為所要求出的根，否則重複(2)以後的步驟。EPS可用收斂判定值來適當的精確度。

26. 2-4 以牛頓法求出方程式的根 • 求xn時，可使用下列的關係式， • 然後以 • 及前面的值xn-1來求出牛頓法的收斂速度比2分法快。

27. class Dr11Panel extends Panel { double f(double x) { return x*x*x-x+1; } double g(double x) { // 微分函數 return 3*x*x-1; } public void paint(Graphics g) { final double EPS=1e-08; // 截斷誤差 final int LIMIT=50; // 截斷次數 double x=-2.0,dx; int k; for (k=1;k<=LIMIT;k++) { dx=x; x=x-f(x)/g(x); if (Math.abs(x-dx)<Math.abs(dx)*EPS) { g.drawString("x="+x,10,20); break; } } if (k>LIMIT) g.drawString("沒有收斂",10,20); } }

28. 2-5 內插 • 例題12 拉格朗基內插設有數組的x,y資料，以拉格朗基內插求出通過這些資料點的內插多項式，並求出資料點以外的點的值。 • 設有的n個點，通過這些點的函數f(x)之求法如下：

29. class Rei12Panel extends Panel { double lagrange(double x[],double y[],int n,double t) { int i,j; double s,p; s=0.0; for (i=0;i<n;i++) { p=y[i]; for (j=0;j<n;j++) { if (i!=j) p=p*(t-x[j])/(x[i]-x[j]); } s=s+p; } return s; } public void paint(Graphics g) { double[] x={0.0,1.0,3.0,6.0,7.0}, y={0.8,3.1,4.5,3.9,2.8}; double t; g.drawString("xy",10,20); for (t=0.0;t<=7.0;t=t+.5) { g.drawString(""+t,10,(int)(t*40+40)); g.drawString(""+lagrange(x,y,5,t),70,(int)(t*40+40)); } } }

30. 2-5 牛頓內插 • 設有 個點，現製作其表如下： • 從此表可求得a0~an-1，以此為係數的n-1次的多項式可由下列式子求出：

31. 2-5 牛頓內插 • 這個計算式稱為牛頓內插多項式。此多項式可使用Horner（練習問題1）方法改寫成： • 求係數a0~an-1時，可使用作業用陣列w[]按下列順序計算。…為陣列w[0]，w[1]…的內容，而w0, w0’, w0”則代表第1次、第2次、第3次…使用相同的陣列w[0]之意。 由於牛頓內插法在一開始會先求出a0~an-1的係數，再以其值為基礎計算出多項式的值後進行內插，因此牛頓內插的計算次數比拉格朗基內插少。

32. class Dr12Panel extends Panel { double newton(double x[],double y[],int n,double t) { int flag=1,i,j; double[] a=new double[100]; // 係數陣列 double[] w=new double[100]; // 作業用 double s; if (flag==1) { // 僅在第1次被呼叫時，才在a[]內求係數 for (i=0;i<n;i++) { w[i]=y[i]; for (j=i-1;j>=0;j--) w[j]=(w[j+1]-w[j])/(x[i]-x[j]); a[i]=w[0]; } flag=-1; } s=a[n-1]; // x=t時的內插 for (i=n-2;i>=0;i--) s=s*(t-x[i])+a[i]; return s; } public void paint(Graphics g) { double[] x={0.0,1.0,3.0,6.0,7.0}, y={0.8,3.1,4.5,3.9,2.8}; double t; g.drawString("xy",10,20); for (t=0.0;t<=7.0;t=t+.5) { g.drawString(""+t,10,(int)(t*40+40)); g.drawString(""+newton(x,y,5,t),70,(int)(t*40+40)); } } }

33. 2-6 多位數計算 • long數與long數的加減法n位數的long數之間的加法與減 法 • 長度無法為long型態、double型態等基本資料型態所處理的的數稱為long數。 • 假設將20位數的數19994444777722229999與01116666333388881111相加，我們不將這2個數存放至1個變數內，而將這2個數由低位數起每4位數分割成1組存放置陣列a[]，b[]內。 然後由a[]與b[]的低位數開始分別加起，再將計算結果存放至c[]內，計算結果如不滿10000，則直接將此計算結果存入c[]內，如為10000以上（超過4位數）則將計算結果減去10000後存入c[]內，在下1個位數相加時再將進位的1加進去。 減法也是以同樣方式進行，計算結果如為負數，則將上10000後存入c[]內，在下1個位數相減時，再將借來的1減去。

34. class Rei13Panel extends Panel { private final int KETA=20, // 位數 N=(KETA-1)/4+1; // 陣列大小 void ladd(int a[],int b[],int c[]) { // long數＋long 數 int i,cy=0; for (i=N-1;i>=0;i--) { c[i]=a[i]+b[i]+cy; if (c[i]<10000) cy=0; else { c[i]=c[i]-10000; cy=1; } } } void lsub(int a[],int b[],int c[]) { // long數－long 數 int i,brrw=0; for (i=N-1;i>=0;i--) { c[i]=a[i]-b[i]-brrw; if (c[i]>=0) brrw=0; else { c[i]=c[i]+10000; brrw=1; } } } void print(Graphics g,int c[],int p) { // 顯示位數 int i; for (i=0;i<N;i++) { String s="0000"+c[i]; // 在前方空白處補上0 g.drawString(s.substring(s.length()-4,s.length()),i*40,p); } } public void paint(Graphics g) { int[] a={1999,4444,7777,2222,9999}, b={ 111,6666,3333,8888,1111}, c=new int[N+2]; ladd(a,b,c);print(g,c,20); lsub(a,b,c);print(g,c,40); } }

35. 2-6 long數與short數的乘除法 • 此處設short數為0～32767的整數。 • long數short數是由低位數起按下列步驟計算： • long數short數是由高位數起按下列步驟計算：

36. class Dr13Panel extends Panel { private final int KETA=20, // 位數 N=(KETA-1)/4+1; // 陣列大小 void lmul(int a[],int b,int c[]) { // long數 x short數 int i,d,cy=0; for (i=N-1;i>=0;i--) { d=a[i]; c[i]=(d*b+cy)%10000; cy=(d*b+cy)/10000; } } void ldiv(int a[],int b,int c[]) { // long數 ÷ short數 int i,d,rem=0; for (i=0;i<N;i++) { d=a[i]; c[i]=(d+rem)/b; rem=((d+rem)%b)*10000; } } void print(Graphics g,int c[],int p) { // 顯示位數 int i; for (i=0;i<N;i++) { String s="0000"+c[i]; // 在前方空白處補上0 g.drawString(s.substring(s.length()-4,s.length()),i*40,p); } } public void paint(Graphics g) { int[] a={ 0,3050,2508,8080,1233}, c=new int[N+2]; lmul(a,101,c);print(g,c,20); ldiv(a,200,c);print(g,c,40); } }

37. 2-7 長的π • 求π的1000位數的正確值 • 根據馬欽(J.Marchin)公式 : • 第n項的表示方式則為 • n若為奇數則為正數，若為偶數，則為負數。 • n用long數計算取得:

38. 2-7 π的多位數計算 … public void paint(Graphics g) { int[] s=new int[L2+2],w=new int[L2+2],v=new int[L2+2],q=new int[L2+2]; int k; for (k=0;k<=L2;k++) s[k]=w[k]=v[k]=q[k]=0; w[0]=16*5;v[0]=4*239; // 馬欽公式 for (k=1;k<=N;k++) { ldiv(w,25,w); ldiv(v,239,v);ldiv(v,239,v); lsub(w,v,q);ldiv(q,2*k-1,q); if ((k%2)!=0) // 判斷為奇數項或偶數項 ladd(s,q,s); else lsub(s,q,s); } printresult(g,s); } } …

39. 2-7 e的1000位數計算 • 練習問題14-1求e的1000位數的正確值 • 根據泰勒展開式，e（自然對數的底）的值可由下列式子求出： • 結果為： 可視情形計算至451項

40. class Dr14_1Panel extends Panel { private final int L=1000, // 所要求出的位數 L1=L/4+1, // 陣列大小 L2=L1+1, // 多取1個 N=451; // 欲計算的項數 ......... void printresult(Graphics g,int c[]) { // 顯示結果 int i,y=0; String s; g.drawString(""+c[0]+".",0,20); // 顯示最高位數 for (i=1;i<L1;i++) { if ((i-1)%15==0) y++; s="0000"+c[i]; g.drawString(s.substring(s.length()-4,s.length()),((i-1)%15)*40,y*20+20); } } public void paint(Graphics g) { int[] s=new int[L2+2],w=new int[L2+2]; int k; for (k=0;k<=L2;k++) s[k]=w[k]=0; s[0]=w[0]=1; for (k=1;k<=N;k++) { ldiv(w,k,w); ladd(s,w,s); } printresult(g,s); } }

41. 2-8 聯立方程式的解法 • 以高斯‧喬丹法解聯立方程式 演算法的內容大致是這樣的： (1)將軸由第1行第1列移至第n行第n列，然後反覆進行下列步驟： (2)以軸係數(akk)除以軸所在行的元素(akk, akk((&+1, …, akn , ak)。軸為1。而軸以前的元素(ak1, ak2+1,…, akk-1)由於已為0，因此不用再除。 (3)軸所在行以外的各行則反覆進行下列步驟： (4)(各行)-(軸所在行) X (係數)。由於軸以前 的列元素已經為0，因此這項步驟可以不用執 行。

42. 以下列的3元聯立方程式為例說明具體的清掃過程：以下列的3元聯立方程式為例說明具體的清掃過程：

43. 2-8　聯立方程式的解法(高斯‧喬丹法) public void paint(Graphics g) { final int N=3; // 原有的數 int i,j,k; double[][] a={{2.0,3.0,1.0,4.0}, // 係數矩陣 {4.0,1.0,-3.0,-2.0}, {-1.0,2.0,2.0,2.0}}; double p,d; for (k=0;k<N;k++) { p=a[k][k]; // 軸係數 for (j=k;j<N+1;j++) // 將軸行除以p a[k][j]=a[k][j]/p; for (i=0;i<N;i++) { // 清掃軸列 if (i!=k) { d=a[i][k]; for (j=k;j<N+1;j++) a[i][j]=a[i][j]-d*a[k][j]; } } } for (k=0;k<N;k++) g.drawString("x"+(k+1)+"="+a[k][N],10,k*20+20); }

44. 2-9 線性規劃法 • 最標準的電腦線性規劃方法為簡單(Simplex)法。其公式為： • 目標函數則設為：

45. 2-9 線性規劃法 加上鬆弛(slack)變數，使其成為等式 再從這項聯立方程式至做出下列的係數矩陣 a11 a12… a1n 1 b1 a21 a22 a2n 1 b2 am1 am2 amn 1 bm -c1 -c2… -cn 0 0 … 0 0

46. 2-9 線性規劃法 • 這種演算法的大致內容如下： (1)從最後1行(目標函數的係數)之中找出最小值(由於此值為負的，因此其絕對值也可說是最大的)所在列y。列的選擇。 (2)如最小值≧0，則結束程式。這項條件意味著最後1行的所有係數為正值，即使再繼續執行清掃運算，目標函數的值也不會再增加。 (3)以(1)所求的y列中的各行元素除以各行右端的元素，再找出最小元素所在行x。行的選擇。 (4)將x行y列設為軸並執行清掃演算。

47. 1 3 1 0 0 18 18 2 3 0 1 0 21 21/2 3 1 0 0 1 21 21/3 ←行的選擇 -2 -1 0 0 0 0 ↑ 列的選擇 0 8/3 1 0 -1/8 11 33/8 0 7/3 0 1 -2/3 7 3 ←行的選擇 1 1/3 0 0 1/3 7 21 0 -1/3 0 0 2/3 14 ↑ 列的選擇 x1 x2 0 0 1 -8/7 3/7 3 0 1 0 3/7 -2/7 3 ←的最適解 1 0 0 -1/7 3/7 6 ←的最適解 0 0 0 1/7 12/21 15 ←z的最適解 z

48. ……… double[][] a={{1.0,2.0,1.0,0.0,0.0,14.0}, // 係數矩陣 {1.0,1.0,0.0,1.0,0.0,8.0}, {3.0,1.0,0.0,0.0,1.0,18.0}, {-2.0,-3.0,0.0,0.0,0.0,0.0}}; double min,p,d; while (true) { min=9999; // 列的選擇 for (k=0;k<N2-1;k++) { if (a[N1-1][k]<min) { min=a[N1-1][k]; y=k; } } if (min>=0) break; min=9999; // 行的選擇 for (k=0;k<N1-1;k++) { p=a[k][N2-1]/a[k][y]; if (a[k][y]>0 && p<min) { min=p; x=k; } } p=a[x][y]; // 軸係數 for (k=0;k<N2;k++) // 以p除以軸行 a[x][k]=a[x][k]/p; for (k=0;k<N1;k++) { // 清掃軸列 if (k!=x) { d=a[k][y]; for (j=0;j<N2;j++) a[k][j]=a[k][j]-d*a[x][j]; } } } for (k=0;k<N3;k++) { flag=-1; for (j=0;j<N1;j++) { if (a[j][k]==1) flag=j; } ……

49. 如要將此方程式與測定資料的距離的2次方和 設為最小，就要決定出係數 這就是最小平方法的考量方向。將右側的聯立方程式解開就可求出最小平方法。而m為適用曲線的階數。 2-10 最小平方法 n組的測定資料 此資料的近似方程式為：

50. 2-10 最小平方法 聯立方程式也能以高斯‧喬丹法來解開，因此我們要製作下列的係數矩陣： 這裡要注意的是，由於係數矩陣的斜線方向為同樣的值，因此無須計算所有的元素。也就是說，各相等的元素內有相同的係數。舉例來說，第2行第0列、第1行第1列、第0行第2列皆有 。