數值計算基礎

1、多項式

1.1 多項式的創建

1.1.1 直接輸入系數向量創建多項式

1.1.2 特殊多項式輸入法

1.1.3 由多項式的根逆推多項式

1.2 多項式的運算

1.2.1 多項式的求值

1.2.2 求多項式的根

1.2.3 多項式的乘除法

1.2.4 多項式的微積分

1.2.5 多項式的部分分式展開

1.2.6 多項式的估值

1.2.7 多項式的擬合

1.2.8 多項式的插值

2、線性方程

2.1 有唯一解線性方程組求法

2.2 有無窮解的線性方程組求法

2.2.1 齊次線性方程組的通解

2.2.2 非齊次線性方程組通解

2.3 自編函式判斷唯一解或者無窮解

3、資料分析

3.1 協方差陣與相關陣

3.1.1 協方差陣

3.1.2 相關陣

3.2 差分和梯度

3.2.1 差分

3.2.2 梯度

1、多項式

1.1 多項式的創建

1.1.1 直接輸入系數向量創建多項式

由于在matlab中多項式的以向量的形式存盤的，直接輸入向量，matlab將降冪自動把向量的元素分配給多項式各項的系數，

創建方法：

在matlab的命令視窗直接輸入多項式的系數矢量，然后利用轉換函式poly2sym將多項式由系數質量形式轉換為符號形式，

例如：輸入系數矢量，創建多項式 $2x^{3}-5x^{2}+3x+7$

>> poly2sym([2 -5 3 7])

ans=

    2x^3-5*x^2+3*x+7

1.1.2 特殊多項式輸入法

創建方法：

matlab提供了poly函式，使用它可以由矩陣的特征多項式創建多項式，使用該方法生成多項式時，其首項的系數必為1.

例如：求矩陣 $\begin{bmatrix} 1 & 2 & 3\\ 4 & 5 & 6\\ 7 & 8 & 0 \end{bmatrix}$ 的特征多項式系數，并轉換為多項式系數，

>> a=[1 2 3; 4 5 6; 7 8 0];
>> p=poly(a)

p=
    1.0000    -6.0000    -72.00000    -27.0000

>> poly2sym(p)

ans=
    x^3-6*x^2-72*x-27

1.1.3 由多項式的根逆推多項式

創建方法：

如果已知某個多項式的根，那么使用poly函式可以很好產生其對應的多項式，

>> roots=[-4 -2+2i -2-2i 5];
>> p=poly(roots)

p=
    1    3    -16    -88    -160
>> disp(poly2sym(p))

x^4+3*x^3-16*x^2-88*x-160

1.2 多項式的運算

1.2.1 多項式的求值

matlab中有指定的函式求多項式的值，呼叫格式如下：

y=polyval(p,x)    % p在x點的值

y=polyval(p,x,[],mu)    % p在x點的值

[y,delta]=polyval(p,x,s)    % 產生誤差估計

y=polyvalm(p,x)    % 其中x是方陣

例如：求解在x=8時多項式(x-1)(x-2) (x-3)(x-4)的值，

   >> p=poly([1 2 3 4]);
   >> polyvalm(p,8)
   ans =
      840

1.2.2 求多項式的根

直接呼叫求根函式roots

那么如何表示系數與根？

》》先把多項式轉化為伴隨矩陣，再求特征值

例如：求解多項式x3-7x2+2x+40的根，

>> r=[1 -7 2 40];

>> p=roots(r);

      -0.2151

       0.4459

       0.7949

       0.2707

1.2.3 多項式的乘除法

c=conv(a,b)%多項式a與b相乘（卷積）

a=deconv(c,b)%多項式相除（解卷積）

>> a=[2 3 4 2 59 8];
   b=[12 2 4 3 5 7 8 9 1];
>> polyval(a,1)%當x=1時，多項式的值
ans =
    78
>>c=conv(a,b)%多項式a與b相乘（卷積）
c=
    24    40    62    50   747   263   315   289   394   508   550   597   131   8
>>deconv(c,b)%多項式相除（解卷積）
d =
    2.0000    3.0000    4.0000    2.0000   59.0000    8.0000
>>roots(a)%當多項式a=0時，x的值
ans =
  -1.8991 + 1.7358i
  -1.8991 - 1.7358i
   1.2172 + 1.7203i
   1.2172 - 1.7203i
  -0.1361 + 0.0000i

例如：計算多項式乘法 $(x^{2}+2x+2)(x^{2}+5x+4)$ ，

>> c=conv([1 2 2],[1 5 4])

   c =

        1     7    16    18     8

計算多項式除法 $(3x^{3}+13x^{2}+6x+8)/(x+4)$ ，

   >> d=deconv([3 13 6 8],[1 4])

   d =

        3     1     2

1.2.4 多項式的微積分

微分運算

k = polyder(p)：求多項式的導函式多項式
k = polyder(a,b)：求多項式a與多項式b乘積的導函式多項式
[q,d] = polyder(b,a)：求多項式b與多項式a相除的導函式，導函式的分子存入q，分母存入d

積分運算

polyint(p,k)：回傳以向量p為系數的多項式積分，積分的常數項為k
polyint(p)：回傳以向量p為系數多項式的積分，積分的常數項為默認值0

例如：計算多項式 $4^{4}-12^{3}-14x^{2}+5x+9$ 的微分和積分，

   >> p=[4 –12 –14 5];

   >> pder=polyder(p);

   >> pders=poly2sym(pder)

   >> pint=polyint(p);

   >> pints=poly2sym(pint)

   pders =

       12*x^2-24*x-14

   pints =

       x^4-4*x^3-7*x^2+5*x

1.2.5 多項式的部分分式展開

matlab中，有理多項式由它們的分子多項式和分母多項式表示，對有理多項式進行運算的兩個函式是residue和polyder，redidue執行部分分式展開的運算

[r,p,k] = residue(b,a)：b,a分別為分子和分母多項式系數的行向量，r為留數行向量
[b,a] = residue(r,p,k)：p為極點行向量，k為直項行向量

例如：對下式進行部分分式展開：

   >> a=[1 3 4 2 7 2];
   >> b=[3 2 5 4 6];
   >> [r,p,k]=residue(b,a)
   r =
      1.1274 + 1.1513i
      1.1274 - 1.1513i
     -0.0232 - 0.0722i
     -0.0232 + 0.0722i
      0.7916          
   p =
     -1.7680 + 1.2673i
     -1.7680 - 1.2673i
      0.4176 + 1.1130i
      0.4176 - 1.1130i
     -0.2991          
   k =
      []

1.2.6 多項式的估值

matlab提供了polyval函式與polyvalm函式用于求多項式p(x)在x=a的取值，輸入可以是標量或矩陣

y = polyval(p,x)：p為多項式的系數向量，x為矩陣，它是按陣列運算規則來求多項式的值
[y,delta] = polyval(p,x,S)：使用可選的結構陣列S產生由polyfit函式輸出的估計引數值；delta是預測未來的觀測估算的誤差標準偏差
y = polyval(p,x,[],mu)或[y,delta] = polyval(p,x,S,mu)：使 $\hat{x} = (x - \mu _{1}) / \mu _{2}$ 替代x， $\mu _{1} = mean(x), \mu _{2} = std(x)$ ，其中心點與坐標值可由polyfit函式計算得出

polyvalm函式的輸入引數只能是N階方陣，這時可以將多項式看作矩陣函式

Y = polyvalm(p,X)：p為多項式的系數向量，X為方陣，其實按矩陣運算規則來求多項式的值，

>> X = pascal(4)
 
X =
 
     1     1     1     1
     1     2     3     4
     1     3     6    10
     1     4    10    20
 
>> p = poly(X)
 
p =
 
    1.0000  -29.0000   72.0000  -29.0000    1.0000
 
>> P = poly2str(p,'x')
 
P =
 
   x^4 - 29 x^3 + 72 x^2 - 29 x + 1
 
>> y = polyval(p,X)
 
y =
 
   1.0e+04 *
 
    0.0016    0.0016    0.0016    0.0016
    0.0016    0.0015   -0.0140   -0.0563
    0.0016   -0.0140   -0.2549   -1.2089
    0.0016   -0.0563   -1.2089   -4.3779
 
>> y = polyvalm(p,X)
 
y =
 
   1.0e-10 *
 
   -0.0003   -0.0036   -0.0052   -0.0143
   -0.0021   -0.0136   -0.0179   -0.0464
   -0.0059   -0.0330   -0.0400   -0.1047
   -0.0130   -0.0639   -0.0750   -0.1962

1.2.7 多項式的擬合

擬合：polyfit()函式，不一定用離散點建立多項式函式關系，但是建立的多項式要與離散點誤差（平方和距離）最小，

p=polyfit(x,y,n)

[p,s]= polyfit(x,y,n)

說明：x,y為資料點，n為多項式階數，回傳p為冪次從高到低的多項式系數向量p，x必須是單調的，矩陣s用于生成預測值的誤差估計，

例如：用多項式擬合法求一個形如 $y=ax^2+bx+c$ 的公式，使它與表中所列資料擬合

>> x=[19 25 31 38 44];
y=[19.0 32.3 49.0 73.3 97.8];
a=polyfit(x,y,2);          %擬合2次函式
x0=19:0.1:44;              %步長為0.1
y0=polyval(a,x0);          %回傳值y0是對應于x0的函式值
plot(x,y,'o',x0,y0,'r')    %畫圖，o表示圓圈，r表示紅色red
legend('擬合前','擬合后')   %給曲線加上說明
xlabel('x');               %給x軸加上說明
ylabel('y'); 
grid on;                   %添加網格線
set(gca,'GridLineStyle',':','GridColor','k','GridAlpha',1);  %將網格線變成虛線
a                          %直接輸出a
 
a =
 
    0.0497    0.0193    0.6882

所以可以得出擬合函式 $y=0.0497x^2+0.0193x+0.6882$

1.2.8 多項式的插值

資料插值可分為多項式（高次）插值、分段（低次）插值、三角插值等，多項式插值包括Lagrange插值、Aitken插值、Newton插值、Hermite插值，但高次插值會出現Runge現象，因此更多使用分段低次樣條插值，

使用最多的為三次樣條插值，

一維多項式插值采用函式interp1( )進行實作
一維快速傅里葉插值通過函式interpft( )來實作，該函式利用傅里葉變換將輸入資料變換到頻域，然后用更多點的傅里葉逆變換，變換回時域，其結果是對資料進行增采樣
二維插值主要用于影像處理和資料的可視化，其基本思想與一維插值相同,對函式進行插值，zi=interp2(x, y, z, xi, yi):通過初始資料x、y和z產生插值函式y= f(x, y)，回傳值zi是(xi, yi)在函式f(x, y)上的值，zi=interp2(x, y, z, xi, yi, method):其中method為可采用的插值方法，二維插值采用的方法只有4種，分別是'nearest'、'linear'. 'spline' 和'cubic',其中線性插值為默認的插值方法，
三次樣條插值可以采用函式spline(),該函式的呼叫格式為:yi=spline(x, y, xi):通過初始資料產生插值函式，然后對資料xi進行插值，回傳值yi=f(xi)，采用這種呼叫方式時,其相當于yi=interp1(x, y, xi, 'spline')，pp=spline(x, y):該函式通過對初始資料x和y產生插值函式，并進行回傳，然后利用函式ppval( )對資料xi進行插值計算,其呼叫方式為yi=ppval(pp, xi)，其中pp為插值函式，

例如：有一正弦衰減資料y=sin(x).*exp(-x/10)，其中x=0:pi/5:4*pi，用三次樣條法進行插值，

    >> x0=0:pi/5:4*pi;  

    >> y0=sin(x0).*exp(-x0/10);

    >> x=0:pi/20:4*pi;

    >> y=spline(x0,y0,x);

    >> plot(x0,y0,'or',x,y,'b')

2、線性方程

2.1 有唯一解線性方程組求法

對于一般的，有唯一解的線性方程組，我們可以轉換成矩陣的形式：A x = b

則可以用矩陣運算求解x，即x=A\b

2.2 有無窮解的線性方程組求法

2.2.1 齊次線性方程組的通解

求解齊次線性方程組基礎解系的函式是null
Z=null(A)表示回傳矩陣A的基礎解系組成的矩陣，Z還滿足 $Z^{T}Z=1$
Z=null(A,‘r’)得出的Z不滿足 $Z^{T}Z=1$ ，但得出的矩陣元素多為整數，顧一般都帶引數r，

2.2.2 非齊次線性方程組通解

非齊次線性方程組在求出基礎決議后還要求一個特解，對于矩陣形式的非齊次線性方程組A x = b 特解x0的求法為x0=pinv(A)*b；其中函式pinv的意思是偽逆矩陣，

例如：解方程組，

   >> a=[2 9 0;3 4 11;2 2 6];

   >> b=[13 6 6]';

   >> x=a\b

   x =

       7.4000

      -0.2000

      -1.4000

求解線性方程組：

由輸出結果可知方程的解為

2.3 自編函式判斷唯一解或者無窮解

撰寫一個函式，給入矩陣A和b，給出方程的解，函式自動判斷是有唯一解還是無窮解，

function varargout = gaussion(varargin)
% gaussion 用高斯消元法解線性方程組
%   A是系數矩陣，b是常熟矩陣，varargin={A,b};如果b為0，則不輸入b
%   varargout=[S flag],S給出結果
%   flag為0無解；1唯一解；2齊次通解；3非齊次通解

A=cell2mat(varargin(1));
Alie=length(A);
Asum=numel(A);
Ahang=Asum/Alie;
if(nargin==2)
    b=cell2mat(varargin(2));
else
    b(Ahang,1)=0;
end
B=A; 
B(:,Alie+1)=b; 


varargout(1)={S};
if(nargout==2)
    varargout(2)={flag};
end
end

Ar=rank(A); Br=rank(B);
B=rref(B);
if (Ar<Br)
    flag=0; S=0;
elseif (Ar==Br && Ar==Alie)
    flag=1; S=B(:,Alie+1);
else
    %將能構成單位矩陣的列號存盤在行向量I中,即存盤了極大線性無關向量的編號
    %將剩余列號存入行向量C中
    for i=1:Ar
        for j=1:Alie
            if(B(i,j)==1)
                I(i)=j;
                break
            end
        end
    end
    C=setdiff(1:Alie,I);
    %由線性代數知識可得基礎解系
    ILim=length(I); CLim=length(C);
    S(Alie,CLim)=0;%初始化S，S行數為A列數，S列數為C的維度
    for i=1:CLim
        S(C(i),i)=-1;
        for j=1:ILim
            S(I(j),i)=B(j,C(i));
        end
    end
    if(nargin==1)
        flag=2;
    else
        flag=3;
        S(Alie,CLim+1)=0;
        for i=1:Ar
            S(I(i),CLim+1)=B(i,Alie+1);
        end
    end
end

3、資料分析

3.1 協方差陣與相關陣

3.1.1 協方差陣

矩陣中的資料按行排列與按列排列求出的協方差矩陣是不同的，這里默認資料是按行排列，即每一行是一個observation(or sample)，那么每一列就是一個隨機變數，

$X_{m\times{n}}=\begin{bmatrix}a_{11} &a_{12} &\cdots &a_{1n} \\ a_{21} &a_{22} &\cdots &a_{2n} \\ \vdots &\vdots & \vdots& \vdots \\a_{m1} &a_{m2} &\cdots &a_{mn} \\ \end{bmatrix}=\begin{bmatrix}c_1 &c_2 &\cdots&c_n \\ \end{bmatrix}$

協方差矩陣：

$covMatrix = \frac{1}{m-1}\begin{bmatrix}cov(c_1,c_1) &cov(c_1,c_2) &\cdots &cov(c_1,c_n) \\ cov(c_2,c_1) &cov(c_2,c_2) &\cdots &cov(c_2,c_n) \\ \vdots &\vdots &\vdots &\vdots \\ cov(c_n,c_1) &cov(c_n,c_2) &\cdots &cov(c_n,c_n) \end{bmatrix}$

clear; clc;

X = [1 2 3;
     3 1 1];
 
Y = X;
 
[rows, cols] = size(X);
 
% get mean of each dimension(each column).
meanMatrix = mean(X);

% X - mean.
X = X - ones(rows, 1) * meanMatrix;

% get the cov matrix.
covMatrix = 1 / (rows - 1) * (X' * X)

% the given 'cov' function
cov(Y)

3.1.2 相關陣

函式 corrcoef

corrcoef(X)：回傳從矩陣X形成的一個相關系數矩陣，若X是一個m*n的矩陣，那么得到的相關系數矩陣A就是一個n*n的對稱矩陣，A中的第i行第j列的元素表示的就是X第i列和第j列的相關系數，

corrcoef(X,Y)：它的作用和corrcoef([X,Y])是一樣的，表示序列x和序列y的相關系數，得到的結果是一個2*2矩陣，其中對角線上的元素分別表示x和y的自相關，非對角線上的元素分別表示x與y的相關系數和y與x的相關系數，兩個是相等的，

corrcoef函式算出來的是皮爾遜相關系數，

corrcoef函式計算相關系數是在matlab提供的cov函式基礎上進行計算的，形成的矩陣是

函式corr

corr(X)輸出的結果和corrcoef是一致的，但是corr可以自己選擇相關系數的型別，matlab提供三種，默認的是皮爾遜相關系數，剩下的兩種是kendall和spearman.

1. X與Y是兩個變數取值所構成的向量

Pearson相關系數：corr(X,Y,'type','Pearson')

Kendall相關系數：corr(X,Y,'type','Kendall')

Spearman相關系數：corr(X,Y,'type','Spearman')

2. X是一個資料矩陣，列為個變數取值

Pearson相關系數：corr(X,'type','Pearson')

Kendall相關系數：corr(X,'type','Kendall')

Spearman相關系數：corr(X,'type','Spearman')

3.2 差分和梯度

3.2.1 差分

計算資料差分的函式diff()

>>  A=[1 2 3;4 5 6;7 8 9; 10 11 12]
 
A =
 
     1     2     3
     4     5     6
     7     8     9
    10    11    12
 
>> B=diff(A,1,1)%按列求一階差分
 
B =
 
     3     3     3
     3     3     3
     3     3     3
 
>> B=diff(A,1,2)%按行求一階差分
 
B =
 
     1     1
     1     1
     1     1
     1     1

3.2.2 梯度

[Fx,Fy]=gradient(F)，其中Fx為其水平方向上的梯度，Fy為其垂直方向上的梯度，Fx的第一列元素為原矩陣第二列與第一列元素之差，Fx的第二列元素為原矩陣第三列與第一列元素之差除以2，以此類推：Fx(i,j)=(F(i,j+1)-F(i,j-1))/2，最后一列則為最后兩列之差，同理，可以得到Fy，

1、如果F是一維矩陣，則FX=gradient(F，H)回傳F的一維數值梯度，H是F中相鄰兩點間的間距，

2、如果F是二維矩陣，回傳F的二維數值梯度，

[FX，FY]=gradient(F,HX,HY) HX，HY引數表示各方向相鄰兩點的距離

3、如果F是三維矩陣，回傳F的三維數值梯度，

[FX，FY，FZ]=gradient(F,HX,HY,HZ) HX，HY，HZ引數表示各方向相鄰兩點的距離，

>> x=[6,9,3,4,0;5,4,1,2,5;6,7,7,8,0;7,8,9,10,0]
x =
     6     9     3     4     0
     5     4     1     2     5
     6     7     7     8     0
     7     8     9    10     0


>> [Fx,Fy]=gradient(x)
Fx =
    3.0000   -1.5000   -2.5000   -1.5000   -4.0000
   -1.0000   -2.0000   -1.0000    2.0000    3.0000
    1.0000    0.5000    0.5000   -3.5000   -8.0000
    1.0000    1.0000    1.0000   -4.5000  -10.0000


Fy =
   -1.0000   -5.0000   -2.0000   -2.0000    5.0000
         0   -1.0000    2.0000    2.0000         0
    1.0000    2.0000    4.0000    4.0000   -2.5000
    1.0000    1.0000    2.0000    2.0000         0

例如：計算運算式的梯度并繪圖

    >> v = -2:0.2:2;
    >> [x,y] = meshgrid(v);
    >> z=10*(x.^3-y.^5).*exp(-x.^2-y.^2);
    >> [px,py] = gradient(z,.2,.2);
    >> contour(x,y,z)
    >> hold on
    >> quiver(x,y,px,py)
    >> hold off

全文共9504個字，碼字總結不易，老鐵們來個三連：點贊、關注、評論

作者：左手の明天

原創不易，轉載請聯系作者并注明出處

敬請期待下一篇：

matlab從無到有系列（四）：符號數學基礎

轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/433373.html

標籤：AI

上一篇：神經網路中的激活函式與損失函式&深入理解推導softmax交叉熵

下一篇：【論文筆記】道路檢測 SNE-RoadSeg

matlab從無到有系列（三）：數值計算基礎

1、多項式

1.1 多項式的創建

1.1.1 直接輸入系數向量創建多項式

1.1.2 特殊多項式輸入法

1.1.3 由多項式的根逆推多項式

1.2 多項式的運算

1.2.1 多項式的求值

1.2.2 求多項式的根

1.2.3 多項式的乘除法

1.2.4 多項式的微積分

1.2.5 多項式的部分分式展開

1.2.6 多項式的估值

1.2.7 多項式的擬合

1.2.8 多項式的插值

2、線性方程

2.1 有唯一解線性方程組求法

2.2 有無窮解的線性方程組求法

2.2.1 齊次線性方程組的通解

2.2.2 非齊次線性方程組通解

2.3 自編函式判斷唯一解或者無窮解

3、資料分析

3.1 協方差陣與相關陣

3.1.1 協方差陣

3.1.2 相關陣

3.2 差分和梯度

3.2.1 差分

3.2.2 梯度