壹、問題
某公交公司1路車過去20個季度內的客流量(單位：百萬)如下表：
季度序號 1 2 3 4 5 6 7 8 9 10
客流量 1.85 2.18 1.6 2.31 1.93 2.35 1.638 2.51 1.92 2.49
季度序號 11 12 13 14 15 16 17 18 19 20
客流量 2.01 2.68 2.05 2.77 2.23 2.95 2.46 2.87 2.35 3.07
試確定客流量與季度序號之間的函數關系，並預測未來八個季度的客流量。
二、實驗目的
掌握最小二乘法的原理並會用於解決實際問題；學會用mathematica中的算符進行曲線擬合。
三、預備知識
1 最小二乘法
在許多實際問題中，往往需要根據實驗測得兩個變量x與y的若幹組實驗數據(x1,y1),…(xn,yn)來建立這兩個變量的函數關系的近似式，這樣得到的函數近似式稱為經驗公式。
通過對實驗數據的處理，能夠判斷x、y大體上滿足某種類型的函數關系y＝f(x,a1,a2,…,as),但是其中s個參數a1,a2,…,as的值需要通過n組實驗數據來確定，通常可以這樣來確定參數：選擇參數a1,a2,…,as，使得f(x,a1,a2,…,as)在x1,x2 …xn處的函數值與實驗數據 y1,y2 …yn 的偏差的平方和為最小，就是使
d＝ (1)
為最小，這種方法稱為最小二乘法。當f(x,a1,a2,…,as)是s個參數的線性函數時，利用求極值與解線性方程組的方法可以解決。
例如，若x、y大體上滿足線性關系即f(x,a,b)＝ax+b ，則
d(a,b)＝ (2)
由多元極值的求法有
(3)
解上述關於a、b的二元壹次方程組得
a＝ ， b＝
從而求得經驗公式y＝ax＋b 。
d＝ 的大小是衡量經驗公式精度的壹種尺度。
線性函數是最簡單最常用的經驗公式，有壹些實際問題，它們的經驗公式可能不是線性函數，我們可以把它化為線性函數來討論，例如y=kemx,兩邊取對數得lny=mx+lnk ,令z=lny,b=lnk,即可化為z＝mx＋b。
2 mathematica中數據擬合算符的用法
在數據處理中常常設法用壹個函數按照某種法則去描述壹組數據，這就是數據擬合。上面介紹的最小二乘法就是壹種最常用的數據擬合方法。mathematica中最基本的數據擬合算符是 fit[ ] ,語法為
fit[數據,擬合函數的基函數列表，變量]
線性函數擬合的基函數為1，x ，n階多項式擬合的基函數是1,x,x2,…xn 。
例 壹冊書的成本費y與印刷的冊數x有關，統計數據如下：
xi(千冊) 1 2 3 4 5 6 7 8 9 10
yi（元） 10.15 5.52 4.08 2.85 2.11 1.62 1.41 1.30 1.21 1.15
試用y=a+ 去擬合上述數據。
mathematica程序及運行結果如下：
data={10.15,5.52,4.08,2.85,2.11,1.62,1.41,1.30,1.21,1.15};
fit[data,{1,1/x},x]
四、實驗內容與要求
1 畫出實驗問題的數據圖，並粗略估計這些數據與什麽類型的函數比較吻合？
2 取經驗公式為線性函數y＝ax＋b 按照最小二乘法的原理用mathematica編程解實驗問題。
3 取經驗公式為y＝ax＋b ＋c sin[ x]+d cos[ x] ,用mathematica中算符fit[]來求解實驗問題，並與內容2的精度比較，對比實際情況，妳能得出什麽？
五、操作提示
1
2 擬合程序及運行結果如下：
預測程序及運行結果如下：
3 程序及運行結果如下：
計算兩種經驗公式的精度可以看出第二種較好，這與客流量呈季節被動變動的實際情況吻合。
怎樣用 mathematica 擬合二元函數？
數據擬合
由壹組已知數據（xk，yk）（k=1，2，…，n），求函數的近似解析式y=f（x），就是數據擬合問題，當然函數還可以是多元的。
Mathematica提供了進行數據擬合的函數：
Fit[data，funs，vars] 對數據data用最小二乘法求函數表funs中各函數的壹個線性組合作為所求的近似解析式，其中vars是自變量或自變量的表。
例如：
Fit[data，{1，x}，x] 求形為y=a+bx的近似函數式。
Fit[data，{1，x，x2}，x] 求形為y=a+bx+cx2的近似函數式。
Fit[data，{1，x，y，x y}，{x，y}] 求形為z=a+bx+cy+dxy的近似函數式。
以上出現的參數data的格式為{{x1，y1，…，f1}，{x2，y2，…，f2}，…}。
函數表中的函數還可以是更復雜的初等函數。
例1 由下面給出的壹組數據進行線性擬合，並繪制擬合曲線。。
xi 19.1 25 30.1 36 40 15.1 50
yi 76.3 77.8 79.25 80.8 82.35 83.9 85.1
解：In[1]：=data={{19.1，76.3}，{25，77.8}，{30.1，79.25}，{36，80.8}，
{40，82.35}，{45.1，83.9}，{50，85.1}}；
f=Fit[data，{1，x}，x]
Out[2]=70.5723+0.291456x
In[3]：= pd=ListPlot[data，DisplayFunction→Identity]；
fd=Plot[f，{x，19，52}，DisplayFunction→Identity]；
Show[pd，fd，DisplayFunction→$DisplayFunction]
圖13-49 線性擬合的示意圖
Out[5]=-Graphics-
說明：上例使用壹次函數得到很理想的擬合，圖形如圖13-49所示。
例2 由下面給出的壹組數據進行二次函數擬合，並繪制擬合曲線。
xi 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
yi 5.1234 5.3057 5.5687 5.9378 6.4337 7.0978 7.9493 9.0253 10.3627
解：In[1]：= data={{0.1，5.1234}，{0.2，5.3057}，{0.3，5.5687}，
{0.4，5.9378}， {0.5，6.4337}，{0.6，7.0978}，
{0.7，7.9493}，{0.8，9.0253}，{0.9，10.3627}}；
f=Fit[data，{1，x，x^2}，x]
Out[2]=5.30661-1.83196x+8.17149x2
In[3]：= pd=ListPlot[data，DisplayFunction→Identity]；
fd=Plot[f，{x，0，1}，DisplayFunction→Identity]；
Show[pd，fd，DisplayFunction→$DisplayFunction]
圖13-50 使用二次函數擬合的示意圖
Out[5]= -Graphics-
以上兩例都是計算方法教材中的習題，利用Mathematica可以輕而易舉地得到答案，並同時畫出圖形以便直觀地了解擬合的質量。
以下是二元擬合。
例3 觀察下面的二元函數擬合。
In[1]：=Flatten[Table[{x，y，1 + 5x –x y}，
{x，0，1，0.2}，{y，0，1，0.2}]，1]
Out[1]={{0，0，1}，{0，0.2，1}，{0，0.4，1}，
{0，0.6，1}，{0，0.8，1}，{0，1.，1}，
{0.2，0，2.}，{0.2，0.2，1.96}，{0.2，0.4，1.92}，
{0.2，0.6，1.88}，{0.2，0.8，1.84}，{0.2，1.，1.8}，
{0.4，0，3.}，{0.4，0.2，2.92}，{0.4，0.4，2.84}，
{0.4，0.6，2.76}，{0.4，0.8，2.68}，{0.4，1.，2.6}，
{0.6，0，4.}，{0.6，0.2，3.88}，{0.6，0.4，3.76}，
{0.6，0.6，3.64}，{0.6，0.8，3.52}，{0.6，1.，3.4}，
{0.8，0，5.}，{0.8，0.2，4.84}，{0.8，0.4，4.68}，
{0.8，0.6，4.52}，{0.8，0.8，4.36}，{0.8，1.，4.2}，
{1.，0，6.}，{1.，0.2，5.8}，{1.，0.4，5.6}，
{1.，0.6，5.4}，{1.，0.8，5.2}，{1.，1.，5.}}
In[2]：=Fit[%，{1，x，y，x y}，{x，y}]
Out[2]=1.+5. x+7.77156×10-16 y -1. x y
In[3]：=Chop[%]
Out[3]= 1.+ 5. x -1. x y
說明：在上例的In[1]中，首先生成二元函數1+5x-xy在0≤x≤1，0≤y≤1時的壹個數據表，然後In[2]由這些數據反過來求二元函數，說明Fit可以求解多元問題。In[3]使用函數Chop去掉系數很小的項，以此消除誤差。
函數Chop的壹般形式為：
Chop[expr，δ] 去掉表達式expr的系數中絕對值小於δ的項，δ的默認值為10-10。
最後這個例子用於說明Fit的第二個參數中可以使用復雜的函數，不限於1，x，x2等基本類型。
例4 觀察下面使用初等函數組合進行的擬合。
In[1]：= ft=Table[N[1+2Exp[-x/3]]，{x，10}]
Out[1]={2.43306，2.02683，1.73576，1.52719，1.37775，
1.27067，1.19394，1.13897，1.09957，1.07135}
In[2]：=Fit[ft，{1，Sin[x]，Exp[-x/3]，Exp[-x]}，x]
Out[2]= 1. -4.44089×10-15e-x +2.e-x/3+2.22045×10-16Sin[x]
In[3]：=Chop[%]
Out[4]=1. +2. e-x/3