function [ ans ] = horner( arr, x )
%输入行向量arr，变量x，求多项式结果p(x)
count = length(arr);
ans = arr(1);
for i=2:count
    ans = ans*x +arr(i);
执行，得到结果和polyval相同：




    

>> polyval(p,3)
ans =
>> horner(p,3)
ans =
二、二次求根
已知二次根式$ax^2+bx+c=0$，求根公式为$(-b\pm\sqrt{b^2-4ac})/(2a)$，公式可化为$(-2c)/(b\pm \sqrt{b^2-4ac})$。
当$|b|≈\sqrt{b^2-4ac}$时，上述式子的分子或分母会变得极小，导致结果造成精度损失，所以，在进行运算时要选择计算相应公式，以减小误差：
  
当b>0，选择（1）的x2和（2）的x1；
  
当b<0，选择（1）的x1和（2）的x2；
同样，我们自己编写一个函数进行计算：
function [ x1 ,x2 ] = FindTheRoot( p )
%输入p = [a b c],求ax^2+bx+c=0的根x1,x2
a = p(1);
b = p(2);
c = p(3);
del = b^2-4*a*c;
x1 = NaN;
x2 = NaN;
if a==0	%一元一次方程
    x1=-c/b;
    x2=-c/b;
    return 
if del < 0	%根不存在
    return 
if b>0 %大于0
    x1 = (-2*c)/(b+sqrt(del));
    x2 = (-b-sqrt(del))/(2*a);
else 	%小于0
    x1 = (-2*c)/(b-sqrt(del));
    x2 = (-b+sqrt(del))/(2*a);
return
测试一下执行结果：
>> [a,b] = FindTheRoot([-1 -4 -2])
   -0.5858
   -3.4142
>> [a,b] = FindTheRoot([1 4 2])
   -0.5858
   -3.4142
三、不动点迭代
不动点指一个实数P，对函数g(x)满足P = g(P)；通过点$p_0$开始，通过$p_{n+1} = g(p_n)$得到一个序列，这一过程称为不动点迭代。
不动点迭代可用于求解方程f(x) = 0，当然，并不是所有的f(x)都可以化为g(x)并求解，需要保证迭代收敛。迭代收敛性的判断不在这里进行说明，下面只给出不动点迭代的Matlab实现：
function [ p, k ,err ,P ] = fixpt( g ,p0 ,tole ,max1 )
%input:
%g：g(x)函数
%p0：第一个迭代点
%tole：允许的误差范围
%max1：最大迭代次数
%output:
%k：迭代次数
%p：所求结果
%err：绝对误差大小
%P：迭代序列
P(1) = p0;
for k = 2 : max1
    p = feval(g,P(k-1));
    P(k) = p;
    err = abs(P(k) - P(k-1));
    relerr = err/(abs(P(k))+eps);
    if err<tole | relerr<tole
        break;
if k == max1
    disp("规定迭代次数内无法得到结果")
然后，我们编写一个函数fun1进行测试：
function y = fun1( x )
y = 1+0.5*sin(x);
这个函数与x = y的交点如图：plot(x,y,x,z,'-.')
利用不动点迭代求解，得到结果如下：
>> [ p, k ,err ,P ] =fixpt('fun1',-1,0.000001,10000)
    1.4987
err =
   1.0668e-06
   -1.0000    0.5793    1.2737    1.4781    
   1.4979    1.4987    1.4987    1.4987
四、二分法求根
对连续函数f(x)，给定一个区间[a,b]，其中f(a)f(b)<0，通过不断二分求函数的根：
function [ c, err ,yc  ] = bisect( f ,a ,b ,tole )
%input:
%f：f(x)函数
%a、b：区间
%tole：允许的误差大小
%output:
%err：绝对误差大小
%yc：f(c)
ya = feval(f,a);
yb = feval(f,b);
if ya*yb>0
    disp("不满足f(a)f(b)<0");
    return
max1 = 1+(log(b-1)-log(tole))/(log(2));	
%终止条件1：最大迭代次数
for k = 1:max1
    c = (a+b)/2;
    yc = feval(f,c);
    if yc== 0 
        err = 0
       return 
    if yb*yc>0
        b = c;
        yb = yc;
        a = c;
        ya = yc;
    if b-a<tole %终止条件2：ab区间大小
        break
c = (a+b)/2;
err = abs(b-a);
yc = feval(f,c);
同样，编写函数fun1 = x^3，，执行：
>> [ c, err ,yc  ] = bisect( 'fun1' ,-2 ,1 ,0.001 )
  -1.2207e-04
err =
   7.3242e-04
yc =
  -1.8190e-12
>> [ c, err ,yc  ] = bisect( 'fun1' ,-2 ,2 ,0.001 )
err =
err =
yc =
五、试值法求根
对连续函数f(x)，给定一个区间[a,b]，其中f(a)f(b)<0，通过取(a,f(a))与(b,f(b))的连线与x轴交点(c,0)，得到序列c，用c拟合f(x)=0的根：
function [ c, err ,yc  ] = regula( f ,a ,b ,tole,epsilon,max1 )
%input:
%f：f(x)函数
%a、b：区间
%tole：c误差大小
%epsilon：f(c)误差大小
%max1：最大迭代次数
%output:
%err：绝对误差大小
%yc：f(c)
ya = feval(f,a);
yb = feval(f,b);
c = NaN;
err = NaN;
yc = NaN;
if ya*yb>0
    disp("不满足f(a)f(b)<0");
    return
for k = 1:max1 %最大迭代次数作为终止条件
    bc = yb*(b-a)/(yb-ya);
    c = b-bc;
    ac = c-a;
    yc = feval(f,c);
    if yc== 0 
        err = 0
       return 
    if yb*yc>0
        b = c;
        yb = yc;
        a = c;
        ya = yc;
    dx = min(ac,bc);%c离a、b的距离最小值
    if dx<tole 
        err = abs(b-a)/2;
        return 
    end %x轴区间作为终止条件
    if abs(yc)<epsilon 
        err = abs(b-a)/2;
        return 
    end %y轴区间作为终止条件
err = abs(b-a)/2;
yc = feval(f,c);
与二分法使用同一函数测例：




    

>> [ c, err ,yc  ] = regula( 'fun1' ,-2 ,1 ,0.001 ,0.001,10000)
    0.1529
err =
    1.0765
yc =
    0.0036
>> [ c, err ,yc  ] = regula( 'fun1' ,-2 ,2 ,0.001 ,0.001,10000)
err =
err =
yc =
六、求解近似根位置
为了粗略估计根在区间中的位置，使用等距采样点将区间分成n个小区间，判断每个区间是否满足根存在的条件：
function R = approot (f,X,epsilon)
%Input
%f：输入函数
%X：等距采样点序列x
%epsilon：y的误差
%Output
%R：近似根序列
n=length(X);
for i=1:n
    Y(i) = feval(f,X(i));
yrange = max(Y)-min(Y);
epsilon2 = yrange*epsilon;
m=0;%近似根数量
X(n+1)=X(n);%防止溢出
Y(n+1)=Y(n);
for k=2:n
if Y(k-1)*Y(k)<=0 %根存在条件1
    m=m+1; 
    R(m)=(X(k-1)+X(k))/2;
s=(Y(k)-Y(k-1))*(Y(k+1)-Y(k));
if (abs(Y(k)) < epsilon2) & (s<=0) %根存在条件2
    m=m+1;
    R(m)=X(k);
测试函数：
>> R = approot ('sin',X,0.00001)
    0.0005    3.1415    6.2835    9.4245
七、牛顿-拉弗森迭代
牛顿法只需要函数f和f的导数，且给定一个起始点，即可开始计算。
function [p0, err, k, y] = newton( f, df, p0, delta, eps, max1 )
%Input
%f：输入函数
%df: f的导数
%p0：初始点
%eps：y的误差范围
%delta：p0的误差范围
%max1：最大迭代次数
%Output
%p0：最终结果
%err：最终误差
%k：迭代次数
%y：f(p0)
for k = 1 : max1
    p1 = p0 - feval(f,p0)/feval(df,p0);
    err = abs(p1-p0);
    p0 = p1;
    y = feval(f,p0);
    if (err<delta)|(abs(y)<eps)
        break;
所用的函数：
function y = f( x )
y = x^3 - 3*x + 2;
function y = df( x )
y = 3*x^2 - 3;
测试（从p0 = 1.2开始）：
>> [p0, err, k, y] = newton('f','df',1.2,0.000000000001,0.0000000001,100)
p0 =
    1.0000
err =
   3.2510e-06
   3.1708e-11
八、割线法
割线法不必计算f的导数，但是需要两个起始点p0、p1，求割线进行迭代。
function [p1, err, k, y] = secant( f, p0, p1, delta, eps, max1 )
%Input
%f：输入函数
%p0、p1：初始点
%eps：y的误差范围
%delta：p0的误差范围
%max1：最大迭代次数
%Output
%p1：最终结果
%err：最终误差
%k：迭代次数
%y：f(p0)
for k = 1 : max1
    y0 =  feval(f,p0);
    y1 =  feval(f,p1);
    p2 = p1 - y1*(p1-p0)/(y1 - y0);
    err = abs(p2-p1);
    p0 = p1;
    p1 = p2;
    y = feval(f,p1);
    if (err<delta)|(abs(y)<eps)
        break;
所用的函数和牛顿法相同，进行测试（从p0 = 1.2开始）：
[p1, err, k, y] = secant('f',1.2,1.21, 0.000000000001,0.0000000001,100)
p1 =
    1.0000
err =
   2.9315e-06
   6.7494e-11