為了證明最近的勤學，貼一些程序。哈哈，實現申明我是matlab菜鳥。不過這次是我一道一道調出來的。
最近學會了直接寫M文件和function函數。
1，最小二乘法（帶權值的）對矩陣進行處理，比較繁瑣，
x=[1 2 3 4 5]';
y=[4 4.5 6 8 8.5]';
w=[2 1 3 1 1]';
A=zeros(3,3);
U=zeros(3,1);
for m=1:5
A(1,1)=A(1,1)+w(m)*sin(x(m))*sin(x(m));
A(1,2)=A(1,2)+w(m)*sin(x(m))*cos(x(m));
A(1,3)=A(1,3)+w(m)*sin(x(m))*exp(x(m));
A(2,2)=A(2,2)+w(m)*cos(x(m))*cos(x(m));
A(2,3)=A(2,3)+w(m)*cos(x(m))*exp(x(m));
A(3,3)=A(3,3)+w(m)*exp(x(m))*exp(x(m));
U(1,1)=U(1,1)+w(m)*sin(x(m))*y(m);
U(2,1)=U(2,1)+w(m)*cos(x(m))*y(m);
U(3,1)=U(3,1)+w(m)*exp(x(m))*y(m);
end
A(2,1)=A(1,2);
A(3,1)=A(1,3);
A(3,2)=A(2,3);
2,復合梯形公式，以及復合拋物線公式
function I=ftrapz(a,b,n) 
h=(b-a)/n;
np1=n+1;
x=a:h:b;
c=ones(np1,1);
c(1,1)=0.5;
c(np1,1)=0.5;
y=exp(-x.^2);
I=h*y*c;
save as ftrapz.m
format long
ftrapz(0,10,1000)

function I=fsimpson(a,b,n) 
h=(b-a)/n;
I1=ftrapz(a,b,n);
a1=a+h*0.5;
b1=b-h*0.5;
x=a1:h:b1;
c=ones(n,1);
fx=exp(-x.^2);
I2=2*h*fx*c;
I=(I1+I2)/3;
save as fsimpson.m
fsimpson(0,10,500)
3：追趕法（雖然有了結果，但是只限于在一道題里面，有function會好很多）
t=cputime;
u=ones(1000,1);
y=ones(1000,1);
l=ones(1000,1);
l(1)=b(1,1);
y(1)=d(1,1)/l(1);
u(1)=c(1,1)/l(1);
for k=2:999
l(k)=b(k,k)-a(k-1,k-1)*u(k-1);
y(k)=(d(k)-a(k-1,k-1)*y(k-1))/l(k);
u(k)=c(k,k)/l(k);
end
k=1000;
l(k)=b(k,k)-a(k-1,k-1)*u(k-1);
y(k)=(d(k)-a(k-1,k-1)*y(k-1))/l(k);
x(1000)=y(1000);
for k=999:-1:1
x(k)=y(k)-u(k)*x(k+1);
4：牛頓法
format long
x(1)=0;
k=2;
x(k)=x(k-1)-(1-x(k-1)-sin(x(k-1)))./(-1-cos(x(k-1)));
while (abs(x(k)-x(k-1))>=1e-3)
k=k+1;
x(k)=x(k-1)-(1-x(k-1)-sin(x(k-1)))./(-1-cos(x(k-1)));
end
disp(x)

5：LU分解
A=[ 1 1 1;-1 3 1;2 -6 1]
b=[6 4 -5]'
[l,u]=lu(A)
y=l^-1*b
x=u^-1*y
y = -5.00000000000000
8.50000000000000
1.50000000000000
x =3
2
1
6：插值
x=[-5:0.1:5]
y=1./(1+x.^2)
plot(x,y,'-b')
hold on
x1=[-5:1:5]
y1=1./(1+x1.^2)
p=polyfit(x1,y1,10)
yy=polyval(p,x1)
plot(x1,yy,'-r')
hold on
x1=[-5:0.5:5]
y1=1./(1+x1.^2)
p=polyfit(x1,y1,20)
yy=polyval(p,x1)
plot(x1,yy,'-y')
hold on
x1=[-5:2.5:5]
y1=1./(1+x1.^2)
p=polyfit(x1,y1,4)
yy=polyval(p,x1)
plot(x1,yy,'-g')
hold on
6：牛頓法解方程組（不會）
[x,y] = solve('x*x*x*x + y*y*y*y = 4','x +y - 2 = 0')
/***************************************************************
* 本算法用最小二乘法依據指定的M個基函數及N個已知數據進行曲線擬和
* 輸入: m--已知數據點的個數M
* f--M維基函數向量
* n--已知數據點的個數N-1
* x--已知數據點第一坐標的N維列向量
* y--已知數據點第二坐標的N維列向量
* a--無用
* 輸出: 函數返回值為曲線擬和的均方誤差
* a為用基函數進行曲線擬和的系數,
* 即a[0]f[0]+a[1]f[1]+...+a[M]f[M].
****************************************************************/
double mini_product(int m,double(*f[M])(double),int n,double x[N],
double y[N],double a[M])
{
double e,ff,b[M][M],c[M][1];
int i,j,k;

for(j=0;j
{
for(k=0;k
{
b[j][k]=0.0;
for(i=0;i
b[j][k]+=(*f[j])(x)*(*f[k])(x);
}
c[j][0]=0.0;
for(i=0;i
c[j][0]+=(*f[j])(x)*y;
}
gaussian_elimination(m,b,1,c); /*求擬和系數*/
for(i=0;i
a=c[0];
e=0.0;
for(i=0;i
{
ff=0.0;
for(j=0;j
ff+=a[j]*(*f[j])(x);
e+=(y-ff)*(y-ff);
}
return(e);
}






/*************************************************************************
* 高斯列主元素消去法求解矩陣方程AX=B,其中A是N*N的矩陣,B是N*M矩陣
* 輸入: n----方陣A的行數
* a----矩陣A
* m----矩陣B的列數
* b----矩陣B
* 輸出: det----矩陣A的行列式值
* a----A消元后的上三角矩陣
* b----矩陣方程的解X
**************************************************************************/
double gaussian_elimination(int n,double a[M][M],int m,double b[M][1])
{
int i,j,k,mk;
double det,mm,f;
det = 1.0;
for(k = 0;k
{
mm=a[k][k];
mk = k;
for(i=k+1;i
{
if(fabs(mm)
{
mm = a[k];
mk = i;
}
}
if(fabs(mm)
return(0);
if(mk!=k) /* 將第K列主元素換行到對角線上*/
{
for(j=k;j
{
f = a[k][j];
a[k][j]=a[mk][j];
a[mk][j]=f;
}
for(j=0;j
{
f = b[k][j];
b[k][j]=b[mk][j];
b[mk][j]=f;
}
det = -det;
}
for(i=k+1;i
{
mm = a[k]/a[k][k];
a[k]=0.0;
for(j=k+1;j
a[j]=a[j]-mm*a[k][j];
for(j=0;j
b[j]=b[j]-mm*b[k][j];
}
det = det*a[k][k];
}
if(fabs(a[k][k])
return 0;
det=det*a[k][k];
for(i=0;i
{
b[n-1]=b[n-1]/a[n-1][n-1];
for(j=n-2;j>=0;j--)
{
for(k=j+1;k
b[j]=b[j]-a[j][k]*b[k];
b[j]=b[j]/a[j][j];
}
}
return(det);
}
§7 MATLAB的應用

　

7.1 MATLAB在數值分析中的應用

插值與擬合是來源于實際、又廣泛應用于實際的兩種重要方法。隨著計算機的不斷發展及計算水平的不斷提高，它們已在國民生產和科學研究等方面扮演著越來越重要的角色。下面對插值中分段線性插值、擬合中的最為重要的最小二乘法擬合加以介紹。

7.1.1 分段線性插值

所謂分段線性插值就是通過插值點用折線段連接起來逼近原曲線，這也是計算機繪制圖形的基本原理。實現分段線性插值不需編制函數程序，MATLAB自身提供了內部函數interp1其主要用法如下：

interp1(x,y,xi) 一維插值

◆ yi=interp1(x,y,xi)

對一組點(x,y) 進行插值，計算插值點xi的函數值。x為節點向量值，y為對應的節點函數值。如果y 為矩陣，則插值對y 的每一列進行，若y 的維數超出x 或 xi 的維數，則返回NaN。

◆ yi=interp1(y,xi)

此格式默認x=1:n ，n為向量y的元素個數值，或等于矩陣y的size(y,1)。

◆ yi=interp1(x,y,xi,’method’)

method用來指定插值的算法。默認為線性算法。其值常用的可以是如下的字符串。

● nearest 線性最近項插值。

● linear 線性插值。

● spline 三次樣條插值。

● cubic 三次插值。

所有的插值方法要求x是單調的。x 也可能并非連續等距的。

正弦曲線的插值示例：

>> x=0:0.1:10;

>> y=sin(x);

>> xi=0:0.25:10;

>> yi=interp1(x,y,xi);

>> plot(x,y,’0’,xi,yi)

則可以得到相應的插值曲線（讀者可自己上機實驗）。

Matlab也能夠完成二維插值的運算，相應的函數為interp2，使用方法與interpl基本相同，只是輸入和輸出的參數為矩陣，對應于二維平面上的數據點，詳細的用法見Matlab聯機幫助。

　

7.1.2 最小二乘法擬合

在科學實驗的統計方法研究中,往往要從一組實驗數據中尋找出自變量x 和因變量y之間的函數關系y=f(x) 。由于觀測數據往往不夠準確，因此并不要求y=f(x)經過所有的點 ,而只要求在給定點上誤差按照某種標準達到最小，通常采用歐氏范數作為誤差量度的標準。這就是所謂的最小二乘法。在MATLAB中實現最小二乘法擬合通常采用polyfit函數進行。

函數polyfit是指用一個多項式函數來對已知數據進行擬合，我們以下列數據為例介紹這個函數的用法：

>> x=0:0.1:1；

>> y=[ -0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2 ]

為了使用polyfit，首先必須指定我們希望以多少階多項式對以上數據進行擬合，如果我們指定一階多項式，結果為線性近似，通常稱為線性回歸。我們選擇二階多項式進行擬合。

>> P= polyfit (x, y, 2)

P=

-9.8108 20.1293 -0.0317

函數返回的是一個多項式系數的行向量，寫成多項式形式為：



為了比較擬合結果，我們繪制兩者的圖形：

>> xi=linspace (0, 1, 100); %繪圖的X-軸數據。

>> Z=polyval (p, xi); %得到多項式在數據點處的值。

當然，我們也可以選擇更高冪次的多項式進行擬合，如10階：

>> p=polyfit (x, y, 10);

>> xi=linspace (0, 1,100);

>> z=ployval (p, xi);

讀者可以上機繪圖進行比較，曲線在數據點附近更加接近數據點的測量值了，但從整體上來說，曲線波動比較大，并不一定適合實際使用的需要，所以在進行高階曲線擬合時，“越高越好”的觀點不一定對的。

　

　

　

　

　

7.2 符號工具箱及其應用

在數學應用中，常常需要做極限、微分、求導數等運算，MATLAB稱這些運算為符號運算。MATLAB的符號運算功能是通過調用符號運算工具箱(Symbolic Math Toolbox)內的工具實現，其內核是借用Maple數學軟件的。MATLAB的符號運算工具箱包含了微積分運算、化簡和代換、解方程等幾個方面的工具，其詳細內容可通過MATLAB系統的聯機幫助查閱，本節僅對它的常用功能做簡單介紹。

7.2.1 符號變量與符號表達式

MATLAB符號運算工具箱處理的對象主要是符號變量與符號表達式。要實現其符號運算，首先需要將處理對象定義為符號變量或符號表達式，其定義格式如下：

格式1： sym (‘變量名’) 或 sym (‘表達式’)

功能： 定義一個符號變量或符號表達式。

例如：

>> sym (‘x’) % 定義變量x為符號變量

>> sym(‘x+1’) % 定義表達式x+1為符號表達式

格式2： syms 變量名1 變量名2 …… 變量名n

功能： 定義變量名1、變量2 ……、變量名 n為符號變量。

例如：

>> syms a b x t % 定義a,b, x,t 均為符號變量

　

7.2.2 微積分運算

1、極限

格式：limit (f, t, a, ‘left’ or ‘right’)

功能：求符號變量t 趨近a 時，函數f 的（左或右）極限。‘left’ 表示求左極限，‘right’ 表示求右極限，省略時表示求一般極限；a省略時變量t 趨近0； t省略時默認變量為x ，若無x則尋找（字母表上）最接近字母x 的變量。

例如：求極限的命令及結果為：

>> syms x t

>> limit ((1+2*t/x)^(3*x) , x, inf )

ans=

exp(6*t)

再如求函數x / |x| ，當時的左極限和右極限，命令及結果為：

>> syms x