matlab.txt

Matlab在化学工程中的应用

为了证明最近的勤学，贴一些程序。哈哈，实现申明我是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