MATLAB学数值分析（一）迭代法解非线性方程

这个学期在学数值分析，课程内容相当于学过的计算方法的升级版，数值分析是一门很有用的学科，可以解决很多工程上实际的问题，学习这门课最好的方法就是把学到的算法自己实现一遍，现在打算开一个新坑，把数值分析每一章学到的算法都用matlab实现一遍。

强力推荐Timothy Sauer的数值分析，他的书没有那么多详细的推导，但是入门举的例子非常恰当，最重要的是他可以帮助你理解算法的本质，知道为什么要这样做，这个系列也会讲一点关于这个算法的本质上的理解

第一篇是关于迭代法求非线性方程

一、二分法

二分法很容易理解，收敛的条件也很简单，满足介值定理即可

function xc=bisect(f,a,b,tol)
if sign(f(a))*sign(f(b)) >= 0
    error('f(a)f(b)<0 not satisfied!') %停止运行
end
fa=f(a);
fb=f(b);
while (b-a)/2>tol
    c=(a+b)/2;
    fc=f(c);
    if fc==0
        break
    end
    if sign(fc)*sign(fa)<0
        b=c;
        fb=fc;
    else
        a=c;
        fa=fc;
    end
end
xc=(a+b)/2;

两个实例
(1)

f=@(x) x+sin(x)-1
bisect(f,1,2,0.001)

解得答案为1.5322

(2)

f=@(x) x+sin(x)-1
bisect(f,1,2,0.001)

解得答案为0.5107

二、不动点法（FPI）

接下来的方法其实都可以归结为迭代法，本质上就是利用所给的方程找出一个形如x=φ(x)的公式，然后进行迭代（可以自行设置迭代精度或者迭代次数来控制求根）

不动点法得到迭代公式的方式十分简单，但是是否收敛以及收敛区间的选取是重点，这里不多提传统的条件，用几何的角度去理解非常形象

例：求x^3+x-1=0的根

有三种不动点迭代的公式

1. x=g(x)=1-x^3
2. x=g(x)=(1-3x)^(1/3)
3. x=g(x)=(1+2x^2)/(1+3x3)

用几何的角度理解就是求y=x与y=g(x)的交点，迭代的过程就是这样一个回旋的过程，这样就不难理解为什么需要g`(x)<L<1
1.x=g(x)=1-x^3
2. x=g(x)=(1-3x)^(1/3)
3.x=g(x)=(1+2x^2)/(1+3*x3)
代码

function [xc,xd]=fpi(g,x0,tol)
x(1)=x0;
x(2)=g(x0);
i=2;
while abs((x(i)-x(i-1))/x(i))>tol
    x(i+1)=g(x(i));
    i=i+1;
end
xc=x(i);
xd=x;

三、牛顿法

牛顿法得到迭代公式的方法可以用泰勒展开来理解

理解了这个之后你可以通过减少忽略的阶数来提高迭代的速度

代码

function f=fun(x)
syms x1 x2;
f1=x1-sin(x1+x2)-1.2;
f2=x2+cos(x1+x2)-0.5;
f=[f1 f2];
end

function df=dfun(x);
f=fun(x);
df=[diff(f,'x1');diff(f,'x2')];
df=conj(df');

function x=newton(x0,tol,N)
con=0;
for i=1:N
    f=eval(subs(fun(x0),{'x1' 'x2'},{x0(1) x0(2)}));
    df=eval(subs(dfun(x0),{'x1' 'x2'},{x0(1) x0(2)}));
    x=x0-f/df;
    for j=1:length(x0)
        il(i,j)=x(j);
    end
    if norm(x-x0)<tol
        con=1;
        break;
    end
    x0=x;
end
%以下是将迭代过程写入txt文档文件名为iteration.txt
fid=fopen('iteration.txt','w');
fprintf(fid,'iteration');
for j=1:length(x0)
    fprintf(fid,'         x%d',j);
end
for j=1:i
    fprintf(fid,'\n%6d     ',j);
    for k=1:length(x0)
        fprintf(fid,' %10.6f',il(j,k));
    end
end
if con==1
    fprintf(fid,'\n计算结果收敛！');
end
if con==0
    fprintf(fid,'\n迭代步数过多可能不收敛！');
end
fclose(fid);

一个实例

newton([1.2 1.2],0.0001,10000)

得到x=1.4336;y=1.4722

迭代过程如下：

iteration x1 x2
1 1.430214 0.792144
2 1.880892 0.992182
3 1.774872 1.773734
4 1.250139 1.864315
5 1.233776 1.506273
6 1.354048 1.183667
7 1.559594 1.114895
8 1.680931 1.423505
9 1.464598 1.726759
10 1.309052 1.657558
11 1.322102 1.432740
12 1.423059 1.272018
13 1.541910 1.312042
14 1.541888 1.517114
15 1.416667 1.630762
16 1.351550 1.553095
17 1.377243 1.414566
18 1.450691 1.347444
19 1.506011 1.410863
20 1.477253 1.529278
21 1.407990 1.564232
22 1.382664 1.499793
23 1.409690 1.420056
24 1.456534 1.401485
25 1.476255 1.456529
26 1.447933 1.518918
27 1.410860 1.521778
28 1.404229 1.475039
29 1.426738 1.433202
30 1.453237 1.435890
31 1.456444 1.474950
32 1.435428 1.504764
33 1.416715 1.496461
34 1.418111 1.465702
35 1.434324 1.446347
36 1.447587 1.455769
37 1.444656 1.480422
38 1.430903 1.492728
39 1.422360 1.482460
40 1.426339 1.463873
41 1.436804 1.456632
42 1.442476 1.466225
43 1.438215 1.480426
44 1.429967 1.484107
45 1.426742 1.475410
46 1.430786 1.465096
47 1.436932 1.463651
48 1.438728 1.471142
49 1.435000 1.478617
50 1.430437 1.478557
51 1.429742 1.472306
52 1.432931 1.467117
53 1.436218 1.467986
54 1.436292 1.473082
55 1.433583 1.476605
56 1.431282 1.475282
57 1.431611 1.471248
58 1.433803 1.468969
59 1.435376 1.470437
60 1.434848 1.473580
61 1.433090 1.474972
62 1.432069 1.473511
63 1.432678 1.471128
64 1.434043 1.470353
65 1.434678 1.471700
66 1.434064 1.473476
67 1.433020 1.473832
68 1.432663 1.472647
69 1.433234 1.471358
70 1.434012 1.471270
71 1.434185 1.472275
72 1.433678 1.473190
73 1.433110 1.473113
74 1.433059 1.472285
75 1.433491 1.471656
76 1.433895 1.471823
77 1.433872 1.472488
78 1.433514 1.472905
79 1.433235 1.472697
80 1.433300 1.472175
81 1.433588 1.471912
82 1.433775 1.472129
83 1.433690 1.472530
84 1.433463 1.472684
85 1.433344 1.472478
86 1.433434 1.472176
87 1.433609 1.472096
88 1.433680 1.472282
89 1.433593 1.472504
90 1.433461 1.472534
91 1.433423 1.472374
92 1.433502 1.472215
93 1.433600 1.472215
计算结果收敛！

四、割线法

割线法则是改进了牛顿法，可以不对目标函数求导便得到迭代公式，将f`(x)替换成f(x(i+1))-f(x(i))/(x(i+1)-x(i))，比较容易理解，不过多讲解，不过这里的代码多加了个一个最高迭代次数的控制，之前的所有方法都可以加一个这个参数。

代码如下：

function [xc,xd]=secant(f,x0,x1,tol,K)
x(1)=x0;
x(2)=x1;
i=2;
while abs((x(i)-x(i-1))/x(i))>tol
    if i>k
    	break;
    end
    x(i+1)=x(i)-f(x(i))*(x(i)-x(i-1))/(f(x(i))-f(x(i-1)));
    i=i+1;
end
xc=x(i);
xd=x;
end

五、练习

最后给大家留两道题自己练习一下

1.steffensen迭代法

steffensen迭代是不动点迭代的一种加速方式，推导过程略，公式如下

尝试用steffensen迭代求方程e^x+10x-2=0的根

2.用迭代法（不动点法）解方程组的根

大家可以关注我的公众号R君的秘密基地回复“数值分析一”就可以得到我的解答啦~