【笔记】函数和方程（1）
目标：理解求方程近似解的二分法、简单迭代法和牛顿迭代法，学会MATLAB内部函数roots、sovle、fsolve、fzero求解方程，并用之解决实际问题。
求方程近似解的简单方法
不存在解析解的方程就需要结合具体方程（函数）的性质，使用作图法或数值法求出近似解。
1、图形法——放大法求根
图形法是分析方程根的性态最简洁的方法。不过，不要总是想得到根的精确值。这些值虽然很粗糙但直观，多少个根，在什么范围，一目了然，并且还可以借助图形局部放大功能，将根定位更加准确一点。
例1 求方程x^5+2x^2+4=0的所有根及大致范围。
解：（1） 首先画出函数f(x)=x^5+2x^2+4的图形，确定方程的实数根的大致范围。
>> clf,ezplot x-x,grid on ,hold on ,ezplot('x^5+2*x^2+4')
回车得出图1-1，由图知，它有一个实根，大致分布在-2与2之间。
（2） 将作图范围不断缩小，用放大法可得到精度越来越高的根的近似值。在MATLAB命令窗口中先后键入：
>> ezplot x-x,grid on ,hold on ,ezplot('x^5+2*x^2+4',[-2,2])
>> ezplot x-x,grid on ,hold on ,ezplot('x^5+2*x^2+4',[-2,-1])
>> ezplot x-x,grid on ,hold on ,ezplot('x^5+2*x^2+4',[-1.6,-1.5])
>> ezplot x-x,grid on ,hold on ,ezplot('x^5+2*x^2+4',[-1.55,-1.54])
可得图1-2。由该图知方程根在-1.544与-1.543之间。
图1
 |
 |
 |
 |
图2
2、数值法
非线性方程f(x)=0求根的方法有区间法和迭代法两大类：二分法、弦位法是区间法；简单迭代法和牛顿法及其变形是迭代法。这里只说二分法、简单迭代法和牛顿法。
1）根的隔离与二分法
根的隔离思想来源于连续函数的零点定理：
若函数f(x)在闭区间[a,b]上连续，且f(a)f(b)<0，则方程f(x)=0在(a,b)区间内至少有一根x*。
二分法是最简单的求根方法，利用连续函数的零点定理，将含根区间逐次减半缩小，取区间的中点构造收敛点列{xn}来逼近根x*。
用该方法求f(x)=0的近似解可以分两步左：
- 确定根的近似位置或大致范围，及确定一个区间[a,b]，是所求根是位于这个区间内的唯一实根。这个区间称为根的隔离区间。（可以通过作图达到）
- 以根的隔离区间为端点作为根的初始近似值，用二分法逐步改进根的近似值的精确度，直至求得满意的近似解。
用二分，理论上区间中点序列{xn}将收敛与根的真值，但收敛速度较慢，所以通常用二分法为其他方法提供初步的近似值。
2）简单迭代法
迭代法的基本原理就是构造一个迭代公式，反复用它得出一个逐次逼近方程根的数列，数列中每一项都是方程根的近似值，只是精度不同。简单迭代法也称逐次迭代法，是非线性方程求根中各类迭代法的基础。
由于对方程做等价变幻根部发生变化，将方程f(x)=0等价变形为x=t(x),构造迭代计算公式x(n+1)=t(x(n))。取定初值x0，算出数列{xn}。如果{xn}收敛与x*，则有
。
这说明，x*就是方程f(x)=0的根。上面x=t(x)称为不动点方程，t(x)称为迭代函数，数列{xn}称为迭代数列。
注：以上t(x)代表
。
3）牛顿迭代法
如果f(x)在[a,b]上具有二阶导数，f(a)f(b)<0，且f’(x)及f’’(x)在[a,b]上保持定号，这时可采用牛顿迭代法求方程f(x)=0在(a,b)内的唯一实根。牛顿迭代法将非线性方程线性化处理为近似方程，然后用近似方程获得求根的迭代方式。具体做法在此不做说明。
牛顿迭代法就是用切线与x轴交点的横坐标近似代替曲线与x轴交点的横坐标。因此牛顿法也称为切线法，是非线性方程求根方法中收敛最快的方法。
收敛定理：
定理 对于方程f(x)=0，若存在区间[a,b]，使得
（1） f(a)f(b)<0，
（2） f’’(x)在[a,b]上连续且保持不变号，
（3） 对任意的x属于[a,b]，f’(x)！=0，
则当迭代初始值x0属于[a,b]，且f’(x)及f’’(x)在[a,b]上保持同号时，牛顿迭代法产生的迭代序列{xn}收敛与方程在区间[a,b]上的唯一解x*。
3、MATLAB中方程求解的基本命令(表1)
| MATLAB命令 | 用途 |
| roots(p) | 求多项式方程的根，其中p是多项式的系数向量 |
| solve(fun) | 求方程fun=0的符号解，如果不能求得精确的符号解，可以计算可变精度的数值解。 |
| solve(fun,var) | 对指定变量var求代数方程fun=0的符号解 |
| fsolve(fun,x0) | 用最小二乘法求非线性方程fun=0在估计值x0附近的近似解 |
| fzero(fun,x0) | 求函数fun在x0附近的零点 |
最小二乘法
1、数据拟合的最小二乘法原理
科学实验或数据统计中，在得不到两组相关数据的精确解析表达式时，就希望找到关系y=F(x)的近似解析式y=F(x)~=f(x)(注：~=代表约等于)；在数据处理中常常设法用一个函数按照某种法则去描述一组数据，这就是数据拟合。这样得到的函数近似式y=f(x)称为经验公式。
根据一组二维数据，及平面上的若干点，要求确定一个一元函数y=f(x)，即曲线，使这些点与曲线总体来说尽量接近，这就是数据拟合称曲线的思想，简称为曲线拟合（Fitting a Curve）。曲线拟合的目的是根据实验获得的数据去建立因变量与自变量之间有效的经验公式，为进一步深入研究提供线索。
给定平面上的点(x1,y1),(x2,y2),…,(xn,yn)，进行曲线拟合有多种方法，其中最小二乘法是最常用的方法。
最小二乘法的原理是：求f(x),使在x1,x2,…,xn处得函数值与实验数据y1,y2,…yn的偏差的平方和为最小。
曲线拟合的实际含义是寻求一个函数y=f(x),使f(x)在某种准则下与所有数据点最为接近，即曲线拟合的最好。最小二乘法就是使所有散点到曲线的距离平方和最小。
已知一组数据，用什么样的曲线拟合最好呢？可以根据散点图进行直观判断，在此基础上，选择几种曲线分别做拟合，然后比较，观察哪条曲线的最小二乘指标最小。
2、MATLAB中曲线拟合的实现
MATLAB软件提供了基本的曲线拟合函数命令：
多项式函数拟合:a=ployfit(xdata,ydata,n)
其中n表示多项式的最高阶数，xdata,ydata为要拟合的数据，要求用数组的方式输入。输出参数a为拟合多项式y=a(1)x^n+…+a(n)x+a(n+1)的系数a=[a(1),…,a(n),a(n+1)]。
多项式在x处得值y可用下面命令计算:
y=polyval(a,x)
例2 求拟合多项式
>> % 给出数据点的x值
>> x=[0.5 1.0 1.5 2.0 2.5 3.0];
>> % 给出数据点的y值
>> y=[1.75 2.45 3.81 4.80 7.00 8.60];
>> % 求出2阶拟合多项式的系数
>> p=polyfit(x,y,2);
>> poly2str(p,'x')
ans =
0.56143 x^2 + 0.82871 x + 1.156
即可用2次多项式f(x)=0.56143 x^2 + 0.82871 x + 1.156拟合已知数据。为了把数据点拟合曲线画在同一张图上，键入
>> % 给出x在0.5~3.0之间的离散值
>> x1=0.5:0.05:3.0;
>> % 求出拟合多项式在x1上的值
>> y1=polyval(p,x1);
>> % 比较拟合曲线效果没得到图2
>> plot(x,y,'*r',x1,y1,'-b')
例3 用二次多项式对悬链线y=1/5ch(x/3)进行拟合，并绘制它们的图形。
解：悬链线方程中，给出x的值计算对应的y是比较困难的，若把函数拟合成二次多项式，计算起来则方便得多。
>> % 将自变量和函数离散化
>> x=-2*pi:0.2:2*pi;
>> y=1/5*cosh(x./3);
>> % 求拟合多项式的系数
>> p=polyfit(x,y,2)
p =
234/15527 -67/226683 301/1639
>> % 求出拟合多项式表达式
>> pk=poly2str(p,'x')
pk =
0.015071 x^2 - 0.00029557 x + 0.18365
>> % 做悬链线函数方程并保留
>> ezplot 1/5*cosh(x/3),hold
Current plot held
>> % 做拟合多项式函数曲线
>> ezplot 0.015071*x^2-0.00029557*x+0.18365
>> grid,legend('悬链线函数线','拟合悬链线二次多项式曲线')
例4 用函数y=a*exp(bx)拟合下列数据
注：线性化拟合的求解结果与非线性拟合有些差异，因为通过对数变换，均方误差的含义已经改变了。
| x | 0.1 | 0.2 | 0.15 | 0 | -0.2 | 0.3 |
| y | 0.95 | 0.84 | 0.86 | 1.06 | 1.50 | 0.72 |
解：若用非线性拟合，记a=c(1),b=c(2),
>> fun=inline('c(1)*exp(c(2)*x)','c','x');
>> x=[0.1 0.1 0.15 0 -0.2 0.3];
>> y=[0.95 0.84 0.86 1.06 1.50 0.72];
>> c=lsqcurvefit(fun,[0,0],x,y)
c =
1.0791 -1.5796
若用线性拟合，y=a*exp(bx)两边取对数得
z=lny=lna+bx
然后通过除法解超定方程组z1=lna+bxi得lna(从而得a)和b。
>> y=[0.95 0.84 0.86 1.06 1.50 0.72];
>> x=[0.1 0.1 0.15 0 -0.2 0.3];
>> z=log([0.95 0.84 0.86 1.06 1.50 0.72]');
>> m=[ones(6,1),x'];
>> c=m\z;
>> a=exp(c(1)),b=c(2)
a =
1.0742
b =
-1.4900
