分享一下数值分析经常遇到的算法，代码有点多；算法原理之类的网上均可以找到，本文只给出对应的代码实现。
1、线性代数的直接接法
%追赶法求解线性方程组Ax=b，其中A是三对角方阵function x=tridiagsolver(A,b)[n,n]=size(A);for i=1:n if(i==1) l(i)=A(i,i); y(i)=b(i)/l(i); else i<n l(i)=A(i,i)-A(i,i-1)*u(i-1); y(i)=(b(i)-A(i,i-1)*y(i-1))/l(i); end if(i<n) u(i)=A(i,i+1)/l(i); endendx(n)=y(n);for j=n-1:-1:1 x(j)=y(j)-x(j+1)*u(j);end
调用程序
clcclearA=[2,-1,0,0;-1,3,-2,0;0,-2,4,-3;0,0,-3,5];b=[6;1;-2;1];x= tridiagsolver(A,b)
2、插值法
2.1 拉格朗日插值法
function yh=lagrange(x,y,xh)n=length(x);m=length(xh);yh=zeros(1,m);for j=1:m; for i=1:n xp=x([1:i-1 i+1:n]); yh(j)=yh(j)+y(i)*prod((xh(j)-xp)./(x(i)-xp)); %注意区分yh和y endend
调用程序
x=[11,12,13];y=[2.3979,2.4849,2.5649];xh=11.75;yh=lagrange(x,y,xh)
2.2 牛顿插值法
function yh=newtonPol(x,y,xh)n=length(x);p(:,1)=x;p(:,2)=y;for j=3:n+1 p(1:n+2-j,j)=diff(p(1:n+3-j,j-1))./(x(j-1:n)-x(1:n+2-j))'; %求差商表 （注意这里有一个 ’ 符号，与差商表不一样的地方）endq=p(1,2:n+1)'; %求牛顿法的系数--取第一行yh=0;m=1;yh=q(1);for i=2:n m=q(i); for j=2:i m=m*(xh-x(j-1)); %求牛顿法中各多项式值(xh-x0)…(xh-xn-1) end yh=yh+m;%求和end
调用程序
x=[11,12,13];y=[2.3979,2.4849,2.5649];xh=11.75;yh= newtonPol(x,y,xh)
3、数值积分与数值微分
3.1 复合梯形公式
function I=ftrapz(f,a,b,n)format long %显示15位双精度h=(b-a)/n;x=linspace(a,b,n+1);y=feval（f,x); I=h*(0.5*y(1)+sum(y(2:n))+0.5*y(n+1));
函数文件
function y=fun1(x)y=exp(-x);
调用程序
t=ftrapz(@fun1,0,1,10)
3.2 复合Simpson公式
function I=fsimpson(fun,a,b,n)h=(b-a)/n;x=linspace(a,b,2*n+1);y=feval（fun,x);I=(h/6)*(y(1)+2*sum(y(3:2:2*n-1))+4*sum(y(2:2:2*n))+y(2*n+1));
函数文件
function y=fun1(x)y=exp(-x);
调用程序
s=fsimpson(@fun1,0,1,10)
4、线性方程组的迭代解法
%高斯-赛德尔迭代法：Function [x,iter]=gs(A,b,tol)D=diag(diag(A));L=D-tril(A);U=D-triu(A);x=zeros(size(b)); %从x=[0;0…]T开始for iter=1:500 x=(D-L)\(b+U*x); %此句换为x=(D)\(b+L*x+U*x);即为Jacobi迭代 error=norm(b-A*x)/norm(b); if(error<tol) break; endend
主函数(调用程序)
A=[2,-1,0;-1,3,-1;0,-1,2];b=[1;8;-5];tol=1e-4;[x,iter]=gs(A,b,tol)
5、非线性方程求根
%牛顿法求非线性方程的根：% 输入：fun--非线性函数；dfun--非线性函数导数；x0--初始值；tol--精度；% 输出：x--非线性方程数值根function [x,iter]=newton(fun,dfun,x0,tol)format longiter=1;x=x0;while iter<500fun,x)/feval（dfun,x); if abs(feval（fun,x))<tol break; end iter=iter+1;end
newton的函数文件
function y=fun3(x)y=x*cos(x)+2;%
newton的导函数文件
function y=dfun3(x)y=cos(x)-x*sin(x);
调用程序
x=newton(@fun3,@dfun3,2,1e-3)
6、矩阵特征值与特征向量的计算
6.1改进乘幂法
function [t,y]=eigIPower(A,v0,ep)[tv,ti]=max(abs(v0));lam0=v0(ti);u0=v0/lam0;err=ep*10; %为第一步循环做准备，此处不考虑0次循环的情况while(err>ep) v1=A*u0; [tv,ti]=max(abs(v1)); lam1=v1(ti); err=abs(lam0-lam1); u0=v1/lam1; lam0=lam1;endt=lam1;y=u0;
调用程序
A=[12,6,-6;6,16,2;-6,2,16];xinit=[1;0.5;-0.5];[t,y]=eigIPower(A,xinit,1e-4)
6.2 反幂法
function [t,y]=eigIPower_inv(A,v0,ep)[tv,ti]=max(abs(v0));lam0=v0(ti);u0=v0/lam0;err=ep*10;while(err>ep) v1=A\u0;
[tv,ti]=max(abs(v1)); lam1=v1(ti); err=abs(1/lam0-1/lam1); %反幂法在误差计算时用的是特征值的倒数 u0=v1/lam1; lam0=lam1;endt=1/lam1;
y=u0;
调用程序
A=[12,6,-6;6,16,2;-6,2,16];xinit=[1;0.5;-0.5];[t,y]=eigIPower_inv (A,xinit,1e-4)
7、常微分方程初边值问题数值解
7.1 标准龙格库塔四阶四段公式
function y=rk4(fun,a,b,y0,n)h=(b-a)/n;y(1)=y0;for k=1:n x=a+(k-1)*h;fun,x,y(k));fun,x+h/2,y(k)+k1/2);fun,x+h/2,y(k)+k2/2);fun,x+h,y(k)+k3); y(k+1)=y(k)+(k1+2*k2+2*k3+k4)/6;end
函数文件
function u=frk4(x,y)u=y-2*x/y;
调用程序
y=rk4(@frk4,0,1,1,10)
7.2 欧拉法
function y=euler(f,a,b,y0,h)n=(b-a)/h;y(1)=y0;for i=1:n x(i)=a+(i-1)*h; y(i+1)=h*feval（f,x(i),y(i));end
函数文件
function u=feuler(x,y)u=x^3-y/x;
调用程序
y=euler(@feuler,1,2,0.4,0.2)