這個學期在學數值分析，課程內容相當於學過的計算方法的升級版，數值分析是一門很有用的學科，可以解決很多工程上實際的問題，學習這門課最好的方法就是把學到的算法自己實現一遍，現在打算開一個新坑，把數值分析每一章學到的算法都用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的根
有三種不動點迭代的公式
- x=g(x)=1-x^3
- x=g(x)=(1-3x)^(1/3)
- x=g(x)=(1+2x^2)/(1+3xxx)
用幾何的角度理解就是求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+3xxx)
代碼
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君的秘密基地回復“數值分析一”就可以得到我的解答啦~
