目錄
- 插值與拟合
- 插值和拟合的差別
- 插值方法
- 分段線段插值
- 拉格朗日插值多項式
- 樣條插值
- 三次樣條插值
- Matlab插值工具箱
- 一維插值函數
- interp1函數
-
- 例題1
- 二維插值
- 例題
- 例題2
- 一維插值函數
- 曲線拟合的線性最小二乘法
- 線性最小二乘法
- 公式推導
- 函數$r_k(x)$的選取
- 最小二乘法的Matlab實作
- 解方程組法
- 例題5.5
- 多項式拟合法
- 線性最小二乘法
- 最小二乘優化
- lsqlin函數
- lsqcurvefit函數
- lsqnonlin函數
- lsqnonneg函數
- Matlab的曲線拟合使用者圖形界面解法
- 曲線拟合與函數逼近
- 曲線拟合
- 函數逼近
圖檔取自知乎使用者yang元祐的回答
插值：函數一定經過原始資料點。
假設f(x)在某區間[a,b]上一系列點上的值
\[y_i=f(x_i),i=0,1,\dots,n。
\]
插值就是用較簡單、滿足一定條件的函數\(\varphi(x)\)去代替\(f(x)\)。插值函數滿足條件
\[\varphi(x_i)=y_i,i=0,1,\dots,n
拟合：用一個函數去近似原函數，不要求過已知資料點，隻要求在某種意義下它在這些點上的總偏差最小。
分線段插值就是将每兩個相鄰的節點用直線連起來，如此形成的一條折線就是就是分段線性插值函數，記作\(I_n(x)\),它滿足\(I_n(x_i)=y_i\),且\(I_n(x)\)在每個小區間\([x_i,x_{i+1}]\)上是線性函數\((i=0,1\dots,n-1)\)。
\(I_n(x)\)可以表示為\(I_n(x)=\sum_{i=0}^n y_il_i(x)\)，其中
\[l_i(x)=
\begin{cases}
\frac{x-x_{i-1}}{x_i-x_{i-1}},&x\in [x_{i-1},x_i],i \neq 0,\\
\frac{x-x_{i+1}}{x_i-x_{i+1}},&x\in [x_i,x_{i+1}],i \neq n,\\
0,&其他
\end{cases}
\(I_n(x)\)有良好的收斂性，即對\(x\in [a,b]\),有
\[\lim _{n \rightarrow \infin}I_n(x)=f(x)
用\(I_n(x)\)計算x點的插值的時候，隻用到x左右的兩個點，計算量與節點個數n無關。但是n越大，分段越多，插值誤差越小。
朗格朗日(Lagrange)插值的基函數為
\[\begin{aligned}
l_i(x)&=\frac{(x-x_0)\cdots(x-x_{i-1})(x-x_{i+1})\cdots(x-x_n)}{(x_i-x_0)\cdots(x_i-x_{i-1})(x_i-x_{i+1})\cdots(x_i-x_n)}\\
&= \prod_{j=0\\j\neq i}^{n} \frac{x-x_j}{x_i -x_j},i=0,1,\cdots,n。
\end{aligned}
\(l_i(x)\)是xn次多項式，滿足
\[l_i(x_j)=
0,&j\neq i,\\
1,& j = i。
拉格朗日插值函數函數
\[L_n(x)=\sum_{i=0}^{n}y_i l_i(x)=\sum_{i=0}^{n} y_i(\prod_{j=0\\j\neq i}^n \frac{x-x_j}{x_i -x_j})
早期工程師制圖時，把富有彈性的細長木條（所謂樣條）用壓鐵固定在樣點上，在其他地方讓它自由彎曲，然後沿木條畫下曲線。成為樣條曲線。繪圖員利用它把一些已知點連結成一條光滑曲線（稱為樣條曲線），并使連接配接點處有連續的曲率。三次樣條插值就是由此抽象出來的。
數學上将具有一定光滑性的分段的分段多項式稱為樣條函數。具體地說，給頂區間[a,b]的一個劃分。
\[\Delta:a=x_0 < x_1 < \dots < x_{n-1} < x_n = b。
- 在每個小區間\([x_i,x_{i=1}](i=0,1,\dots,n-1)\)上是S(x)是m次多項式。
- S(x)在[a,b]上具有m-1階連續函數。
則稱S(x)為關于劃分\(\Delta\)的m次樣條函數，其圖形為m次樣條函數。
已知函數\(y=f(x)\)在區間[a,b]上的n+1個節點
的值\(y_i=f(x_i)(i0,1,\dots,n)\)，求插值函數S(x)，使得
- \(S(x_i)=y_i(i=0,1,\dots,n)\)
- 在每個小區間\([x_i,x_{i+1}](i=0,1,\dots,n-1)\)上S(x)是三次多項式，記為\(S_i(x)\)
- \(S_i(x)\)在[a,b]上二階連續可微。
由條件2,我們記
\[S(x)=\left \{ S_i(x),x\in[x_i,x_{i+1}],i=0,1,\dots,n-1 \right \}\\
S_i(x)=a_i x^3+b_i x^2+c_i x + d_i，
\(a_i,b_i,c_i,d_i\)為待定系數，共4n個
由條件3中得二階連續可微，有
\[\begin{cases}
S_i(x_{i+1})=S_{i+1}(x_{i+1}),\\
S_i^{'}(x_{i+1})=S_{i+1}^{'}(x_{i+1}),i=0,1,\dots,n-2。\\
S_i^{''}(x_{i+1})=S_{i+1}^{''}(x_{i+1}),\\
由上面的式子共确定4n-2方程，為确定S(x)的4n個參數，常用的确定三次樣條函數邊界條件有3種類型：
-
\(S'(a)=y_0',S(b)'=y_n'\)，由這種邊界條件建立的樣條插值函數稱為f(x)的完備三次樣條插值函數。
特别的，\(y_0'=y_n'=0\)時，樣條曲線呈水準狀态。
如果\(f'(x)\)不知道，我們可以使\(S'(x)\)與\(f'(x)\)在端點處近似相等。這時以\(x_0,x_1,x_2,x_3\)為節點作一個三次Newton插值多項式\(N_a(x)\)。同理，以\(x_n,x_{n-1},x_{n-2},x_{n-3}\)為節點作一個三次Newton插值多項式\(N_b(x)\),要求
\[S'(a)=N'_a(a),S'(b)=N'_b(b)。
由這種邊界條件建立的三次樣條稱為\(f(x)\)的Lagrange三次樣條插值函數。
- \(S''(a)=y''_0,S''(b)=y''_n\)。特别地，\(y''_0=y''_n=0\)時，稱為自然邊界條件
- \(S'(a+0)=S'(b-0),S''(a+0)=S''(b-0)\)。此條件稱為周期條件。
y = interp1(x0,y0,x,'method')
% method 為插值方法，預設為線性插值，其值可為
% 'nearest' 最近項插值
% 'linear' 線性插值
% 'spline' 立方樣條插值
% 'cubic' 立方插值
所有的插值方法要求x0是單調的。
當x0為等距時可以使用快速插值法，使用快速插值法的格式為
*nearest
,
*linear
*spline
*cubic
以下為matlab的官方說明
vq = interp1(x,v,xq)
vq = interp1(x,v,xq,method)
vq = interp1(x,v,xq,method,extrapolation)
vq = interp1(v,xq)
vq = interp1(v,xq,method)
vq = interp1(v,xq,method,extrapolation)
pp = interp1(x,v,method,'pp')
說明
vq = interp1(x,v,xq)
使用線性插值傳回一維函數在特定查詢點的插入值。向量 x 包含樣本點，v 包含對應值 v(x)。向量 xq 包含查詢點的坐标。
如果您有多個在同一點坐标采樣的資料集，則可以将 v 以數組的形式進行傳遞。數組 v 的每一列都包含一組不同的一維樣本值。
vq = interp1(x,v,xq,method)
指定備選插值方法：'linear'、'nearest'、'next'、'previous'、'pchip'、'cubic'、'v5cubic'、'makima' 'spline'。預設方法為 'linear'。
vq = interp1(x,v,xq,method,extrapolation)
用于指定外插政策，來計算落在 x 域範圍外的點。如果希望使用 method 算法進行外插，可将 extrapolation 設定為 'extrap'。您也可以指定一個标量值，這種情況下，interp1 将為所有落在 x 域範圍外的點傳回該标量值。
vq = interp1(v,xq)
傳回插入的值，并假定一個樣本點坐标預設集。預設點是從 1 到 n 的數字序列，其中 n 取決于 v 的形狀：
- 當 v 是向量時，預設點是 1:length(v)。
- 當 v 是數組時，預設點是 1:size(v,1)。
如果您不在意點之間的絕對距離，則可使用此文法。
vq = interp1(v,xq,method)
指定備選插值方法中的任意一種，并使用預設樣本點。
vq = interp1(v,xq,method,extrapolation)
指定外插政策，并使用預設樣本點。
pp = interp1(x,v,method,'pp')
使用 method 算法傳回分段多項式形式的 v(x)。
interp1官方文檔
Matlab種資料點稱為斷點。如果三次樣條插值沒有邊界條件，最常用的方法，就是采用非扭結(not - a -kont)條件。這個條件強迫第1個和第2個三次多項式的三階導數相等。對最後一個和倒數第2個多項式也做相同的處理。
% matlab中三次樣條插值有以下函數
y = interp1(x0,y0,x,'spline');
y = spline(x0,y0,x);
pp = csape(x0,y0,conds);
pp = csape(x0,y0,conds,valconds);y=fnval(pp,x);
% x0, y0是已知資料點；x是插值點，y是插值點的函數值
對于三次樣條插值，推薦使用函數
csape
，
csape
的傳回值是pp形式，要獲得插值點的函數值，必須調用函數
fnval
,即為
pp = csape(x0,y0,conds,valconds);y=fnval(pp,x);
pp = csape(x0, y0);% 預設邊界條件，Lagrange邊界條件
pp = csape(x0, y0, conds, valconds);
% valconds 設定邊界的二階導數值為[0,0]
% conds指定插值的邊界條件,其值可為
% 'complete' 邊界我為一階導數，一階導數的值在valconds參數中給出，若忽略valconds參數，按預設情況處理
% 'not - a - knot' 非扭結條件
% 'periodic' 周期條件
% 'second' 邊界為二階導數，二階導數的值在valconds參數中給出，若忽略valconds參數，按預設情況處理
對于特殊條件，可以通過conds的一個\(1 \times 2\)矩陣來表示，conds的取值為0,1,2
例如，conds=[2,1]的意思為，左邊界是二階導數，右邊界是一階導數。對應的值由valconds給出。
csape官方文檔
如下
| t | 0.15 | 0.16 | 0.17 | 0.18 |
|---|---|---|---|---|
| v(t) | 3.5 | 1.5 | 2.5 | 2.8 |
用三次樣方插值求位移\(S=\int_{0.15}^{0.18}v(t)dt\)
clc;
clear;
x0=[0.15,0.16,0.17,0.18];
y0=[3.5,1.5,2.5,2.8];
pp=csape(x0,y0); % 預設的邊界條件，Lagrange邊界條件
format long g
xinshu = pp.coefs; % 顯示每個區間上三次多項式的系數
s = quadl(@(t)ppval(pp,t),0.15,0.18); % 求積分
format % 恢複短小數的顯示格式
若節點是二維的，插值函數就是二進制函數，即曲面。
Matlab中計算二維插值的指令，如：
z = interp2(x0,y0,z0,x,y,'method')
如果是三次樣條插值，可以使用指令
pp = csape({x0,y0},z0,conds,valconds),z=fnval(pp,{x,y})
interp2官方文檔
高程資料點
| y \ x | 100 | 200 | 300 | 400 | 500 |
|---|---|---|---|---|---|
| 636 | 697 | 624 | 478 | 450 | |
| 698 | 712 | 630 | 420 | ||
| 680 | 674 | 598 | 412 | ||
| 662 | 626 | 552 | 334 | 310 |
Q:找出最高點和該點的高程。
clc;
clear;
x = 100:100:500;
y = 100:100:400;
z = [636,697,624,478,450;
698,712,630,478,420;
680,674,598,412,400;
662,626,552,334,310];
pp = csape({x,y},z');
xi = 100:10:500;
yi = 100:10:400;
cz = fnval(pp,{xi,yi});
[i,j]= find(cz==max(max(cz)));
% 要用兩層max，因為max(cz)為y=180時，和x=100:10:500的一系列值，max(max(cz))才是z的最大值。
x = xi(i);
y = yi(j);
zmax = cz(i,j);
>> [x,y]
ans =
170 180
>> zmax
zmax =
720.6252
海底水深資料
| x | 129 | 140 | 103.5 | 88 | 185.5 | 195 | 105 | 157.5 | 107.5 | 77 | 81 | 162 | 117.5 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| y | 7.5 | 141.5 | 23 | 147 | 22.5 | 137.5 | 85.5 | -6.5 | -81 | 3 | 56.5 | -66.5 | 84 | -33.5 |
| z | 4 | 8 | 6 | 9 |
Q:繪制海底曲面的圖形
clc;
clear;
x = [129,140,103.5,88,185.5,195,105,157.5,107.5,77,81,162,162,117.5];
y = [7.5,141.5,23,147,22.5,137.5,85.5,-6.5,-81,3,56.5,-66.5,84,-33.5];
z = -[4,8,6,8,6,8,8,9,9,8,8,9,4,9];
xmm = minmax(x);
ymm = minmax(y);
xi = xmm(1):xmm(2);
yi = ymm(1):ymm(2);
zi1 = griddata(x,y,z,xi,yi','cubic');% 立方插值
zi2 = griddata(x,y,z,xi,yi','nearest'); % 最近點插值
% 立方插值和最近點插值的混合插值的初始值
zi = zi1;
zi(isnan(zi1))=zi2(isnan(zi1));% 把立方插值中不确定值換成最近點插值的結果
subplot(1,2,1),plot(x,y,'*');
subplot(1,2,2),mesh(xi,yi,zi);% 繪制三維圖形
注：Matlab插值時外插值是不确定的，這裡使用了混合插值，把不确定的插值換成了最近點插值的結果。
\[f(x)=a_1r_1(x)+a_2r_2(x)+\dots+a_mr_m(x)
\(r_k(x)\)為事先選好的x一組線性無關的函數;\(a_k\)為待定系數\((k=1,2,\dots,m;m<n)\)。
定義：最小二乘法就是\(y_i(k=1,2,\dots,n)\)與\(f(x_i)\)的距離\(\delta_i\)的平方和最小，是以稱為最小二乘法
\[J(a_1,\dots,a_m)=\sum_{i=1}^n\delta_i^2=\sum_{i=1}^{n}[f(x_i)-y_i]^2
利用取得極值的必要條件\(\frac{\partial J}{\partial a_j}=0\),得到關于\(a_1,\dots,a_m\)的線性方程組,即分别對每一個a求偏導。
\[\sum_{i=1}^n r_j(x_i)\left[ \sum_{k=1}^{m}a_kr_k(x_i)-y_i \right]=0,j=1,\dots,m,
即，
\[\sum_{i=1}^n a_k\left[ \sum_{k=1}^{m}r_j(x_i)r_k(x_i)\right]= \sum_{k=1}^{m}r_j(x_i)y_i,j=1,\dots,m。
記
\[R=
\left[
\begin{matrix}
r_1(x1) & \dots & r_m(x_1)\\
\vdots & \vdots & \vdots\\
r_1(x_n) & \cdots & r_m(x_n)\\
\end{matrix}
\right]_{n\times m}\\
A=[a_1,\cdots,a_m]^T,Y=[y_1,\cdots,y_n]^T,
方程組可以表示為
\[R^T RA=R^TY。
當\(\left \{ r_1(x),\cdots,r_m(x) \right \}\)線性無關時，R列滿秩，\(R^TR\)可逆，于是
\[A=(R^TR)^{-1}R^TY
函數\(r_k(x)\)的選取
常用的曲線有
- 直線\(y=a_ix+a_2\)
- 多項式\(y=a_1x^m+\cdots+a_mx+a_{m+1}\)（一般m=2,3，不宜太高）
- 雙曲線(一支)\(y=\frac{a_1}{x}+a_2\)
- 指數曲線\(y=a_1e^{a_2x}\),
- 對于指數曲線，拟合前需作變量替換，化為對a1,a2的線性函數
選取時，可在直覺判斷的基礎上，選幾種曲線分别拟合，然後比較，選擇最小二乘名額J最小的一個。
記為
\[J(a_1,\dots,a_m)=\Vert RA-Y \Vert_2^2
Matlab中線性最小二乘的标準型為
\[\min_A \Vert RA-Y \Vert_2^2
指令為
A = R\Y
Q:用最小二乘法求一個形如\(y=a+bx^2\)的經驗公式，使其與下列資料表拟合
| 19 | 25 | 31 | 38 | 44 | |
|---|---|---|---|---|---|
| 32.3 | 49.0 | 73.3 | 97.8 |
clc;
clear;
x = [19,25,31,38,44]';
y = [19.0,32.3,49.0,73.3,97.8]';
r = [ones(5,1),x.^2];
ab=r\y;
x0=19:0.1:44;
y0=ab(1)+ab(2)*x0.^2;
plot(x,y,'o',x0,y0,'r')
如果取\(\left \{ r_1(x),\cdots,r_{m+1} \right \}=\left \{ 1,x,\cdots,x^m \right \}\)，即用m次多項式來拟合給定資料。
Matlab指令
a = polyfit(x0,y0,m)
其中，x0,y0為要拟合的資料；m為對項式的次數。輸出參數a為拟合多項式\(y=a(1)x^m+\cdots+a(m)x+a(m+1)\)的系數向量\(a=[a(1,),\cdots,a(m),a(m+1)]\)
求多項式在x處的值y可用以下指令
y = polyval(a,x)
我們用多項式拟合來拟合上面的例題
clc;
clear;
x = [19,25,31,38,44]';
y = [19.0,32.3,49.0,73.3,97.8]';
a = polyfit(x,y,2);
xi = 19:0.1:44;
yi = polyval(a,xi);
plot(x,y,'o',xi,yi,'r')
如果我們比較一下兩者的差別
clc;
clear;
x = [19,25,31,38,44]';
y = [19.0,32.3,49.0,73.3,97.8]';
a = polyfit(x,y,2);
xi = 19:0.1:44;
yi = polyval(a,xi);
r = [ones(5,1),x.^2];
ab=r\y;
x0=19:0.1:44;
y0=ab(1)+ab(2)*x0.^2;
plot(x0,y0,xi,yi)
legend('最小二乘','多項式拟合')
我們看到其實兩者差别不大的，如果我們看一看系數
| \(a(n)\) | a(1) | a(2) | a(3) |
|---|---|---|---|
| 多項式拟合 \(y_1\) | 0.0497 | 0.0193 | 0.6881 |
| 最小二乘法 \(y_2\) | 0.0500 | 0.9725 |
\[y_1=0.0497x^2+0.0193x+0.6881\\
y_2=0.9725x^2+0.05
在無限制優化問題中，有些情形，比如目标函數由若幹個函數的平方和構成，這類函數一般可以寫成
\[F(x)=\sum_{i=1}^mf_i^2(x),x\in R^n,
式中，\(x=[x1,\cdots,x_n]^T\),一般假設\(m\geq n\)。
把極小化這類函數的問題
\[\min F(x)=\sum_{i=1}^mf_i^2(x)
稱為最小二乘優化問題。
在Matlab優化工具箱中，有
lsqlin, lsqcurvefit, lsqnonlin, isqnonneg等函數
求解
\[\min _x \frac{1}{2} \Vert C \cdot x -d \Vert_2^2\\
s.t.
A \cdot x \leq b,\\
Aeq \cdot x = beq,\\
lb \leq x \leq ub,
式中，C, Aeq, A為矩陣；d, b, beq, lb, ub, x為向量
Matlab中的函數為
x = lsqlin(C,d,A,b)
x = lsqlin(C,d,A,b,Aeq,beq,lb,ub)
x = lsqlin(C,d,A,b,Aeq,beq,lb,ub,x0,options)
x = lsqlin(problem)
[x,resnorm,residual,exitflag,output,lambda] = lsqlin(___)
%lsqlin指令求解例5.5
clc;
clear;
x = [19,25,31,38,44]';
y = [19.0,32.3,49.0,73.3,97.8]';
r = [ones(5,1),x.^2];
ab=lsqlin(r,y);
x0=19:0.1:44;
y0=ab(1)+ab(2)*x0.^2;
plot(x,y,'o',x0,y0,'r')
計算結果是一樣的
lsqlin官方文檔
給定輸入輸出數列xdata,ydata，求參量x，使得
\[\min _x \Vert F(x,xdata)-ydata \Vert_2^2 = \sum_i(F(x,xdata_i)-ydata_i)^2
x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)
% fun為定義函數F(x,xdata)的M檔案
注：非線性拟合時，每一次的運作結果可能都是不同的。
Q:用最小二乘法拟合\(y=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^2}{2\sigma^2}}\),其中未知參數為\(\sigma,\mu\)
clc;
clear;
x0 = -10:0.01:10;
y0 = normpdf(x0,0,1);
mf=@(cs,xdata)1/sqrt(2*pi)/cs(2)*exp(-(xdata-cs(1)).^2/cs(2)^2/2);
% yc = mf([2,1],1);% 測試匿名函數
cs = lsqcurvefit(mf,rand(2,1),x0,y0);% 拟合參數的初始值時任意取的
% 計算出來的估計值 cs(1)=0,cs(2)=1
lsqcurvefit官方文檔
已知函數向量\(F(x)=[f_1(x),\cdots,f_k(x)]^T\)，使x使得
\[\min _x \Vert F(x) \Vert_2^2
x = lsqnonlin(fun,x0,lb,ub,options)
% fun為定義向量函數F(x)的M檔案
lsqnonlin官方文檔
求解非負的x，使得，
\[\min _x \Vert Cx-d \Vert_2^2
x = lsqnonneq(C,d,options)
lsqnonneq官方文檔
Matlab工具箱提供了指令
cftool
，該指令給出了一維資料拟合的互動式環境。
執行步驟：
- 把資料導入到工作空間
- 運作
cftool
- 選擇适當的模型進行拟合
- 生成一些相關的統計量
曲線拟合是已知一組離散資料\(\left \{ (x_i,y_i),i=1,\cdots,n \right \}\)，選擇一個較簡單的的函數f(x)（如多項式），在一定的準則（如最小二乘法準則）下，最接近這些資料。
如果已知一個較為複雜的連續函數\(f(x),x\in [a,b]\)，要求選擇一個較簡單的函數f(x)，在一定的準則下最接近f(x)，就是所謂的函數逼近
與最小二乘準則相對應，函數逼近常采用的一種準則是最小平方逼近
\[J=\int_a^b [f(x)-y(x)]^2dx
達到最小。與曲線拟合一樣，選一組函數\(\left \{ r_k(x),k=1,\cdots,m \right \}\)構造函數f(x)，即令
\[f(x)=a_1r_1(x)+\cdots+a_mr_m(x),
帶入上式中，求\(a_1,\cdots,a_m\)使J達到最小。利用極值必要條件可得
\[\left[
(r1,r_1) & \cdots & (r_1,r_m)\\
\vdots & \vdots & \vdots \\
(r_m,r_1) &\vdots & (r_m,r_m)\\
\right]
a_1\\
\vdots\\
a_m
=
(y,r_1)\\
\vdots \\
(y,r_m)
\right],
這裡\((g,h)=\int_a^b g(x)h(x)dx\)，當方程組的系數矩陣非奇異時，有唯一解。
最簡單的是使用多項式逼近，\(r_1(x)=1,r_2(x)=x,r_3(x)=x^2,\cdots\)。并且如果能使\(\int_a^b r_i(x)r_j(x)dx=0,i \neq j\)，方程組的系數矩陣将是對角陣，計算大大簡化，滿足這種性質的多項式稱為正交多項式。
勒讓德(Legendre)多項式是在[-1,1]區間上的正交多項式，它的表示式為
\[P_0(x)=1,P_k(x)=\frac{1}{2^kk!}\frac{d^k}{dx^k}(x^2-1)^k,k=1,2,\cdots
可以證明
\[\int_{-1}^1 P_i(x)P_j(x)dx=
0,&i \neq j,\\
\frac{2}{2i+1},&i=j,
\\
P_{k+1}(x)=\frac{2k+1}{k+1}xP_k(x)-\frac{k}{k+1}P_{k-1}(x),k=1,2,\cdots。
常用的正交多項式還有第一類切比雪夫(Chebyshev)多項式
\[T_n(x)=\cos(narccosx),x\in [-1,1],n=0,1,2,\cdots
和拉蓋爾(Laguerre)多項式
\[L_n(x)=e^x \frac{d^n}{dx^n}(x^ne^{-x}),x\in [0,+\infin),n=0,1,2,\cdots
clc;
clear;
syms x
base=[1,x^2,x^3];
y1 = base.'*base;
y2 = cos(x) *base.';
r1 = int(y1,-pi/2,pi/2);
r2 = int(y2,-pi/2,pi/2);
a = r1\r2;
xishu1=double(a); % 符号資料轉化成數值型資料
xishu2=vpa(a,6); % 把符号資料轉化為保留6位有效數字的符号資料