数值分析部分算法Octave代码
Chapter 2 Solutions of Equations int One VariableAlgorithm 2.1 The Bisection Method (二分法)CodeAlgorithm 2.2 Fixed-Point Iteration (不动点迭代)CodeAlgorithm 2.3 Newton's Method (Newton迭代法)CodeAlgorithm 2.4 Secant (正割法)CodeAlgorithm 2.5 False Position (试位法)CodeChapter 3 Lagrange Interpolating PolynomialsAlgorithm 3.1 Neville's Iterated Interpolation (Neville迭代插值法)Algorithm 3.3 Neville's Interpolation (Neville插值)CodeAlgorithm 3.5 Natural Cubic Spline (自然三次样条)Chapter 4 Numerical DifferentiationAlgorithm 4.1 Composite Integration: Simpson's Rule Algorithm (复合Simpson法则)CodeAlgorithm 4.2 RombergCodeChapter 5 Elementary Theory of Initial-Value ProblemsAlgorithm 5.1 Euler's MethodCode参考书：Numerical Analysis (9th Edition) R L Burder & J D Faires
代码没用Matlab写，用的是Online Octave免费的，google搜一下Octave就能找到。
Chapter 2 Solutions of Equations int One Variable
一元方程求解
Algorithm 2.1 The Bisection Method (二分法)
Code
clc; clear;close all;% LEI YUXIN Blog:/Cynthia018function val = func(x)val = x - 2*sin(x);endtol = 1.e-5;a = 1.0;b = 2.0;nmax = 100;% Initializationiteration = 0;error = 1.0;% iteration begins herefprintf('%-10s %-10s %-10s %-10s %-10s %-8s %s\n','iteration','a','b','p','fa','fp','error');while (iteration <= nmax && error >= tol)iteration++; % Generate and save iteratresp = a + (b - a) / 2;z(iteration) = p;fa = func(a);fp = func(p);error = abs(a-p)/(abs(p) + eps);fprintf('%5d %10.5f %10.5f %10.5f %10.6f %10.7f %10.6f\n',iteration,a,b,p,fa,fp,error);if (error < tol)p_final = p;elseif (fa*fp < 0) % root is between a and pb = p;else %root is between p and ba = p;endend
Algorithm 2.2 Fixed-Point Iteration (不动点迭代)
Code
clc;clear;close all;% LEI YUXIN Blog:/Cynthia018function val=g(x)val=(2-exp(x)+x^2)/3;endtol = 1.e-5;nmax = 100;%initializationp = [0]; % initialized p0count = 1; %iteration = count-1abs_error = [1.0];while (count <= nmax && abs_error(count) >= tol)p(count+1) = g(p(count));abs_error(count+1) = abs(p(count+1)-p(count));count ++; % Generate and save countendif ( count-1 < nmax)%print tableprintf(" n |\tpn\t| abs.error |\n");printf("----------------------------------\n");for i=1:countprintf(" %d | %.7f\t| %-10.6f|\n",i-1,p(i),abs_error(i));endprintf("The number of iterations is %d", count-1); elseprintf("Exceeded maximum iterations");end
Algorithm 2.3 Newton’s Method (Newton迭代法)
Code
clc;clear;close all;% LEI YUXIN Blog:/Cynthia018function val=f(x)val=cos(x)-x;endfunction dval=df(x)dval=-sin(x)-1;% dval=diff(g,x);endtol = 1.e-5;nmax = 100;%initializationp = [1]; % initialized p0count = 1; %iteration = count-1abs_error = [1.0];while (count <= nmax && abs_error(count) >= tol)p(count+1) = (p(count) - f(p(count))/df(p(count)));abs_error(count+1) = abs(p(count+1)-p(count));count ++; % Generate and save countendif ( count-1 < nmax)% print tableprintf(" n | p(n-1)\t| pn\t\t| pn-p(n-1)\t| abs.error |\n");printf("----------------------------------\n");for i=2:countprintf(" %d | %.7f\t| %.7f\t| %.7f\t| %-10.6f|\n",i-1,p(i-1),p(i),p(i)-p(i-1),abs_error(i));endprintf("The number of iterations is %d", count-1); elseprintf("Exceeded maximum iterations");end
Algorithm 2.4 Secant (正割法)
Code
clc;clear;close all;% LEI YUXIN Blog:/Cynthia018function val=f(x)val=cos(x)-x;endtol = 1.e-5;nmax = 100;abs_error(2) = [1.0];count = 2;p = [0.5];q0 = f(p(1));p(2) = pi/4;q1 = f(p(2));while (count < nmax && abs_error(count) >= tol)p(count+1) = p(count) - q1*(p(count)-p(count-1))/(q1-q0);abs_error(count+1) = abs(p(count+1)-p(count));count++;q0 = q1;q1 = f(p(count));endif ( count-1 < nmax)%print tableprintf(" n |p(n-2)\t|p(n-1)\t\t|pn\t\t| abs.error |\n");printf("----------------------------------\n");for n=3:countprintf(" %d | %.7f\t| %.7f\t| %.7f\t| %-10.6f|\n",n-1,p(n-2),p(n-1),p(n),abs_error(n));endprintf("The number of iterations is %d", count-1); elseprintf("Exceeded maximum iterations");end
Algorithm 2.5 False Position (试位法)
Code
clc;clear;close all;% LEI YUXIN Blog:/Cynthia018function val=f(x)val=cos(x)-x;endtol = 1.e-5;nmax = 100;abs_error(2) = [1.0];count = 2;p = [0.5];q0 = f(p(1));p(2) = pi/4;q1 = f(p(2));p0 = p(1);p1 = p(2);while (count < nmax && abs_error(count) >= tol)p(count+1) = p1 - q1*(p1-p0)/(q1-q0);abs_error(count+1) = abs(p(count+1)-p1);count++;q = f(p(count));if (q*q1 < 0)p0 = p1;q0 = q1;endp1 = p(count);q1 = q;endif ( count-1 < nmax)%print tableprintf(" n |p(n-2)\t|p(n-1)\t\t|pn\t\t| abs.error |\n");printf("----------------------------------\n");for n=3:countprintf(" %d | %.7f\t| %.7f\t| %.7f\t| %-10.6f|\n",n-1,p(n-2),p(n-1),p(n),abs_error(n));endprintf("The number of iterations is %d", count-1); elseprintf("Exceeded maximum iterations");end
Chapter 3 Lagrange Interpolating Polynomials
Algorithm 3.1 Neville’s Iterated Interpolation (Neville迭代插值法)
Algorithm 3.3 Neville’s Interpolation (Neville插值)
Code
clc;clear;close all;% LEI YUXIN Blog:/Cynthia018%initializationx=[0.32, 0.35]; %x0,x1,...,xnf=[0.3146, 0.3429]; %f(x)f2=[0.9492, 0.9394]; %f'(x)z=[];Q=[];%initialize z and Q arraysfor i=1:size(x,2)%duplicate xz(2*i-1)=x(i);z(2*i)=x(i);%duplicate f(x)Q(2*i-1,1)=f(i);Q(2*i,1)=f(i);end%starting with j=2 as first column is already initializedfor j=2:size(z,2) %size of array z in 2nd dimension (number of columns)for i=j:size(z,2)%if i is an even number, and j is 2 (2nd column)%is where z0=z1, hence use derivative of f(x)if (mod(i,2)==0 && j==2)Q(i,j)=f2(i/2);elseQ(i,j)=(Q(i,j-1)-Q(i-1,j-1))/(z(i)-z(i-j+1));endendend%display Qprintf("Array of Q:\n");for i=1:size(Q,1)for j=1:size(Q,2)printf("%10f ",Q(i,j));endprintf("\n");end
Algorithm 3.5 Natural Cubic Spline (自然三次样条)
Chapter 4 Numerical Differentiation
Algorithm 4.1 Composite Integration: Simpson’s Rule Algorithm (复合Simpson法则)
Code
clc; clear;close all;% LEI YUXIN Blog:/Cynthia018function val=f(x)val=exp(x);enda=0;b=4;%h=2;h=1;n=(b-a)/h;XI0=f(a)+f(b);XI1=0;XI2=0;for i=1:n-1X=a+i*h;if(mod(i,2)==0)XI2=XI2+f(X);elseXI1=XI1+f(X);endendXI=h*(XI0+2*XI2+4*XI1)/3;printf("If ");hprintf("XI=%f",XI);
Algorithm 4.2 Romberg
Code
clc; clear;close all;% LEI YUXIN Blog:/Cynthia018function val=f(x)val=sin(x);enda=0;b=pi;h=b-a;n=6;R=[];%R(1,1)=h/2*(f(a)+f(b));R(1,1)=0;printf("%.8f\n",R(1,1));%step3for i=2:n%step4sum=0;for k=1:2^(i-2)sum=sum+f(a+(k-0.5)*h);endR(2,1)=1/2*(R(1,1)+h*sum);printf("%.8f\t",R(2,1));%step5for j=2:iR(2,j)=R(2,j-1)+(R(2,j-1)-R(1,j-1))/(4^(j-1)-1);printf("%.8f\t",R(2,j));endprintf("\n");%step7h/=2;for j=1:iR(1,j)=R(2,j);endend
Chapter 5 Elementary Theory of Initial-Value Problems
Algorithm 5.1 Euler’s Method
Code
clc; clear;close all;% LEI YUXIN Blog:/Cynthia018a=0;b=2;h=0.5;N=(b-a)/h;t=a;alpha=0.5;w=alpha;printf("t\t\tw\n");printf("%f\t%f\n",t,w);for i=1:Nw=w+h*(w-((i-1)*h)^2+1);t=a+i*h;printf("%f\t%f\n",t,w);end
