[自high] 數值分析筆記 更新至2022.12.26

Smallcat

目次：
數學模型(mathematical modeling) 2F
數值方法(numerical methods) 2F
matlab相關 2F
根與優化(root and optimization) 3F
線性系統(linear systems) 4F
特徵值(eigenvalue) 6F
曲線擬合(curve fitting) 7F
傅立葉分析(fourier analysis) 8F
多項式差值(polynomial interpolation) 9F
樣條及分段插值(splines and piecewise interpolation) 10F

Smallcat

數學模型(mathematical modeling)

簡單數學模型：
Dependent variable = f(independent variable, parameters, forcing functions)
parameters取決於系統本身的性質。
forcing functions是指外在影響的作用力。

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
matlab基本操作

向量使用方法：
>> a = [1 2 3 4 5]
a = 1 2 3 4 5
>> b = [2;4;6;8;10] or a = [2 4 6 8 10]'
b = 2
      4
      6
      8
      10
>> c = [1 2 3; 4 5 6; 7 8 9]
c = 1 2 3
     4 5 6
     7 8 9
>> b(4)
ans = 8
>> c(2,3)
ans = 6

>> E = zeros(2,3)
E = 0 0 0
      0 0 0
>> u = ones(1,3)
u = 1 1 1

數列使用方法：
>> t = 1:5
t = 1 2 3 4 5
>> t = 1:0.5:3
t = 1.0 1.5 2.0 2.5 3.0
>> t = 10:-1:5
t = 10 9 8 7 6 5

linspace：
產生平均分布的row向量，格式為linspace(x1,x2,n)。
>> linspace(0,1,6)
ans = 0.0 0.2 0.4 0.6 0.8 1.0

四捨五入(round)、最近上數整數(ceil)、最近下數整數(floor)。

graphics(圖形)：這個我真的不會用，求教學
plot(t, v, '0')。
subplot(m, n, p)，將不同資料分散開來顯示。

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
M-files

函式寫法：
function outvar = funcname(arglist)
statements
outver = value;
end

global X
x = value;

關係operators：
==相等、~=不相等、<、>

邏輯運算：
~x x&y x||y



動畫(animation)：
for j=1:n
plot commands
M(j) = getframe;
end
movie(M)
movie(m,n,fps)，m存放frames，n重複次數預設1。

匿名函數(anonymous func)：
格式：fhandle = @(arglist) expression
>> f1 = @(x,y) x^2+y^2
>> f1(3,4)
ans = 25

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
捨入截斷誤差(round off and truncation errors)

true value = approximation + error

true fractional relative error = (true value - approximation)/true value
or(因為我們不知道真實值)
approximation fractional relative error = (present approximation - previous approximation)/present approximation

捨入誤差：
數位電腦無法精確表示某些數量，因此導致誤差。

截斷誤差：
近似的過程中把高次項忽略導致的誤差。

Smallcat

根與優化(root and optimization)

最基礎就是作圖法找交點。
透過變號與否偵測函數圖形是否通過y=0。

包圍法(bracketing methods)：

二分法(bisection)，隨著疊代次數縮小搜尋範圍。
中點和兩點之中，必定恰有一點和中點符號相異，以此作為基準逐漸縮小範圍。


false position


開放法(open methods)：不需要包圍根，但有可能發散

fixed-point iteration
若存在實數使得g(a)=a，則稱a為為函數g的fixed-point。
令Xi+1 = g(Xi)，以此疊代。

Newton-Raphson method
根據微分近似定義延伸的方法。
Xi+1 = Xi - f(Xi)/f'(Xi)，以此疊代。

正割法(secant methods)
若是不方便微分的情況，將牛頓法的微分項改成割線近似微分。
雖然用到兩個x的估計值，但不牽涉正負號轉換，所以不算包圍法。
Xi+1 = Xi - f(Xi)*(Xi-1 - Xi)/(f(Xi-1)-f(Xi))，以此疊代。

brent's method(不太會，需要求救)
上述幾種算法和插值法的組合。

黃金分割搜尋(搜尋最小值)：
已知黃金比例φ為φ^2-φ-1=0的正數解，約等於1.618。
X1 = Xlow + d
X2 = Xup - d
d = (φ-1)(Xup-Xlow)
縮小範圍的邏輯為，若f(X1)<f(X2)，x2成為新的Xlow，x1為新的x2，同時為該次的Xopt；
若f(X2)<f(X1)，x1成為新的Xup，x2為新的x1，同時為該次的Xopt。
判別是否停止疊代的error為(2-φ)*|(Xup-Xlow)/Xopt|。

Smallcat

高斯消去(Gaussian Elimination)

因此第一步我們要先建立好新的係數，用下方的雙重迴圈完成。

再來就是按照上面x的公式解出來。


克拉瑪法則(Cramer's rule)

LU分解(LU factorization)：
L：下三角矩陣，U：上三角矩陣。
Ax = b => LUx = b => Ld = b(求d) => Ux = d(求x)
利用高斯消去法可以做出L和U。

先求d，d的公式如下。

最後根據d求出x，公式如下。


科列斯基分解(Cholesky factorization)
將對稱矩陣A寫成[U]'[U]的形式，其值會符合下方形式。

matlab用來找U的function叫做chol。

Smallcat

疊代法

高斯-賽德爾疊代(Gauss–Seidel method)
雅可比法(Jacobi Method)

Smallcat

Eigenvalue (λ)
(A-λI)X = 0
matlab內建輸入A求特徵方程式的功能：
A = [10 -5; -5 10];
p = poly(A) => p = 1 -20 75
意思是λ^2 - 20λ + 75。
代入roots()渴求出λ的值。

power method
一種疊代逼近法，可決定最大或最重要的eigenvalue，稍微調整後也能找出最小eigenvalue。
找尋的過程中能夠順便找到對應eigenvector。
以下例題先假設X=[1 1 1]，以做出來的結果歸一化後當作新的X。

以前後X的變化量作為Error估計，若是能收斂的話最終error會越來越小。

matlab有直接算出eigenvalue和eigenvector的function「eig」。
e = eig(A)，回傳A矩陣所有eigenvalue。
[V D] = eig(A)，V為eigenvector，D為相對應的eigenvalue。

Smallcat

曲線擬合(curve fitting)

線性回歸(Linear regression)
令離均差平方和為St。
標準差Sy = (St/(n-1))^0.5，matlab的函數為std()。
變異係數(Coefficient of variation) = Sy/y平均

r = rand(m, n) 生成m-by-n的隨機矩陣(0~1)。
r = randn(m, n) 生成m-by-n的隨機矩陣，會常態分佈且平均值為0、標準差為1。
利用線性轉換，r = mn + s * rand(m, n)能生出平均值為mn、標準差為s的隨機矩陣。

線性最小平方回歸(least squares method)
設回歸線 y = ao + a1*x +e，e為回歸線與各資料點的誤差，故e = y - a0 - a1*x，目標為令e最小，回歸線最合適。
Sr = sum(e^2)，令誤差的平方和為Sr。(誤差有正有負，平方後便沒有這個問題)
求Sr極小值，可微分求根，最後推導出a0和a1為下圖公式。

因同時考慮到x和y，自由度-2，所以最小平方法的標準差為(Sr/(n-2))^0.5。

令r^2 = (St - Sr)/St，r為決定係數。
期望中Sr越小越好，所以當Sr=0、r^2=1，代表該最小平方回歸線完美符合。
當St = Sr、r^2=0代表最小平方沒有比較有幫助。
下面是電腦應用的r公式。


多項式回歸(polynomial regression)
前面提過Sr，這邊推廣成最高次二次，即Sr = sum((yi - ao - a1*xi - a2xi^2)^2)。
將Sr分別對三個a做偏微分，重新整理後會得到下式。

若是知道a以外的數據，就能利用矩陣除法(或是說反矩陣乘法)求得三個a，進而求出Sr並算出標準差，要注意每多一個次方就減少自由度1。
使用polyfit(x,y,n)這個func也能算出a，n為最高次數。

一般型態線性最小平方(general linear least squares)
令y = a0z0 + a1z1 + ... + amzm + e，可寫成{y} = [Z]{a} + {e}，其中z是下方的矩陣。

和前面同樣邏輯，Sr可寫成下方的樣子。

因為[Z]a = y，所以[Z]T*[Z]a = [Z]T*y，a = [Z]T*[Z] \ [Z]T*y。

非線性回歸(nonlinear regression)
設f(a0,a1) = sum((y - ~~~)^2)。
使用matlab內建函數，[x, fval] = fminsearch(fun, x0, options, p1...)。
x為讓fun最小的參數，fval為fun的最小值，x0為事前猜測的參數，作為後續推演的初始值使用。
option沒有要特別選的話，填入[]，最後的p為代入fun的數據，如果option為[]，系統會自己找比較適合代入的位置。

Smallcat

弦波最小平方
y = A0 + A1cos(wt) + B1sin(wt) + e，此式的cos和sin項可以合併成一個sin或cos。
Sr = sum(y - (A0 + A1cos(wt) + B1sin(wt))^2)
使用以前的微分解法，因為週期函數的因素，有些項會被消掉，最後會得到以下結果。


連續傅立葉數列(continuous fourier series)
f(t) = a0 + Σ[akcos(kwt) + bksin(kwt)]


快速傅立葉轉換(FFT)
內建函數F = fft(f, n)，F是一個包含DFT的向量，f式包含信號的向量，n是選擇性參數，象徵使用者想使用幾個點的FFT。

Smallcat

條件數(Condition number)
matlab內建cond(A)計算出Ax=b的問題有多容易計算，條件數越大代表b的誤差會越大。

牛頓多項式差值(Newton Interpolating Polynomial)
最簡單的形式，兩個資料點有一個直線，兩點的斜率相等，以此寫出下式並轉換。

二次差值的話可以寫成f2(x) = b1 + b2(x-x1) + b3(x-x1)(x-x2)。
因此各係數為：
b1 = f(x1)
b2 =  [f(x2) - f(x1)]/(x2-x1)。

以此方式可以推廣到n次，其過程可以用迴圈方便電腦計算，用疊代的方式一層一層計算(如下)


下面有課本範例的牛頓多項式插值程式：

1. function yint = Newtint(x,y,xx)
   
2. % Newtint: Newton interpolating polynomial
   
3. % yint = Newtint(x,y,xx): Uses an (n - 1)-order Newton
   
4. % interpolating polynomial based on n data points (x, y)
   
5. % to determine a value of the dependent variable (yint)
   
6. % at a given value of the independent variable, xx.
   
7. % input:
   
8. % x = independent variable
   
9. % y = dependent variable
   
10. % xx = value of independent variable at which
    
11. % interpolation is calculated
    
12. % output:
    
13. % yint = interpolated value of dependent variable
    
14. 
    
15. % compute the finite divided differences in the form of a
    
16. % difference table
    
17. n = length(x);
    
18. if length(y)~=n, error('x and y must be same length'); end
    
19. b = zeros(n,n);
    
20. % assign dependent variables to the first column of b.
    
21. b(: 1) = y(:); % the (:) ensures that y is a column vector.
    
22. for j = 2:n
    
23. 　　for i = 1:n-j+1
    
24. 　　　　b(i,j)=(b(i+1,j-1)-b(i,j-1))/(x(i+j-1)-x(i));
    
25. 　　end
    
26. end
    
27. 
    
28. % use the finite divided differences to interpolate
    
29. xt = 1;
    
30. yint = b(1,1);
    
31. for j = 1:n-1
    
32. 　　xt = xt*(xx-x(j));
    
33. 　　yint = yint+b(1,j+1)*xt;
    
34. end

複製代碼

拉格朗日插值法(Lagrange interpolating polynomial)
將f(x)拆解成∑Li(x)f(xi)，其中Li(x)=Π(x-xj)(xi-xj)。
下面有課本範例的拉格朗日多項式插值程式：

1. function yint = Lagrange(x,y,xx)
   
2. % Lagrange: Lagrange interpolating polynomial
   
3. % yint = Lagrange(x,y,xx): Uses an (n - 1)-order
   
4. % Lagrange interpolating polynomial based on n data points
   
5. % to determine a value of the dependent variable (yint) at
   
6. % a given value of the independent variable, xx.
   
7. % input:
   
8. % x = independent variable
   
9. % y = dependent variable
   
10. % xx = value of independent variable at which the
    
11. % interpolation is calculated
    
12. % output:
    
13. % yint = interpolated value of dependent variable
    
14. 
    
15. n = length(x);
    
16. if length(y)~=n, error('x and y must be same length'); end
    
17. S = 0;
    
18. for i = 1:n
    
19. 　　product = y(i);
    
20. 　　for j = 1:n
    
21. 　　　　if i ~= j
    
22. 　　　　　　product = product*(xx-x(j))/(x(i)-x(j));
    
23. 　　　　end
    
24. 　　end
    
25. 　　s = s+product;
    
26. end
    
27. yint = s;

複製代碼

Smallcat

樣條函數
將高階多項式拆成許多低階片段多項式的連結，這樣的結果被稱為樣條函數(spline function)。
假設有n個資料點，就會有n-1個區間，每個區間有一個獨立的樣條函數si(x)。
現在我們令這些樣條函數都是線性的，所以可以設si(x) = ai + bi(x-xi)。
其中ai = fi，bi = (fi+1 - fi)/(xi+1 - xi)
因此si(x) = fi + [(fi+1 - fi)/(xi+1 - xi)]*(x - xi)
當然這些樣條也可以是二次或是更多次，只要方便就好(?)

尋找表(Lookup Table)
尋找表是用簡單的查詢操作替換執行時計算的陣列，在使用前會先建立幾筆資料，要使用的時候找到接近的答案後插值。

二次樣條(quadratic splines)
二次的樣條能確保有連續的一次、二次微分，這些是最頻繁使用的。
二次樣條si(x) = ai + bi(x - xi) + ci(x - xi)^2
假設現在有n個資料點，n-1個樣條，就有3(n-1)個未知數，因此需要3(n-1)個方程式來解。
第一步，連續條件，函數要經過每個點，si(x)代入fi時，b、c兩項為0，所以ai=fi。
第二步，相鄰的兩個樣條會有共同點(x=xi+1)，即fi + bi(x - xi) + ci(x - xi)^2 =fi+1 + bi+1(x - xi+1) + ci+1(x - xi+1)^2
令hi=xi+1 - xi，可化簡為fi+bihi+cihi^2=fi+1。
第三步，共同點微分後要相等，即-bi - 2ci(x - xi) = -bi+1 - 2ci+1(x - xi+1)。
可化簡為bi + 2cihi = bi+1。
上面兩步可以解出大部分的b和c。
第四步，二次樣條的二次微分等於2ci，令第一點的二次微分為0，即ci=0。

三次樣條(cubic splines)
同理當然也有三次樣條，si(x) = ai + bi(x - xi) + ci(x - xi)^2 + di(x - xi)^3。
a的部分一樣，ai = fi。
bcd的部分類似，fi + bihi + cihi^2 +dihi^3 = fi+1。
分別微分一次和兩次後和原式三條式子可解出大部分bcd，應用上一章的f[xi,xj]定義可化簡為下式。

這個式子可以解出i從2~n-2的bcd。
假設起點和終點的二次微分為0，分別解出c1 = -3d1h1和cn-1 + 3dn-1hn-1 = 0。(其中cn特別設為0)
以上的條件可以看成下方矩陣的乘法。


matlab內建片段插值函數
yy = spline(x, y, xx)，x和y向量儲存被插值的點值，xx向量為輸入的變數，yy向量儲存依據xx資料進行樣條插值的結果。
若是y比x多兩個值，首值會被視為首點的一次微分值，尾值會被視為末點的一次微分值。

yi = interp1(x, y, xi, 'method')，x和y向量儲存被插值的點值，xi向量為輸入的變數，yi向量儲存依據xi資料進行樣條插值的結果。
method選擇想要的插值方法，有'linear'(預設)、'nearest'、'next'、'previous'、'pchip'、'cubic'、'v5cubic'、'makima' 或 'spline'。