常用数值算法Matlab源码(PDE见另贴)

[原创]常用数值算法Matlab源码（PDE见另贴）【更新完成】数值计算方法中的一些常用算法的Matlab源码，这些程序都是原创，传上来仅供大家参考，不足之处请大家指正，切勿用作商业用途……说明：这些程序都是脚本函数，不可直接运行，需要创建函数m文件，保存时文件名必须与函数名相同，懂一点儿Matlab的朋友应该知道。每个程序的说明里面都附了测试例子这些程序在Vista操作系统下，使用MatlabR2008b(7.7版本)编写1、Newdon迭代法求解非线性方程function[xkt]=NewdonToEquation(f,df,x0,eps)%牛顿迭代法解线性方程%[xkt]=NewdonToEquation(f,df,x0,eps)%x:近似解%k:迭代次数%t:运算时间%f:原函数，定义为内联函数%Df:函数的倒数，定义为内联函数%x0:初始值%eps:误差限%%应用举例：%f=inline('x^3+4*x^2-10');%df=inline('3*x^2+8*x');%x=NewdonToEquation(f,df,1,0.5e-6)%[xk]=NewdonToEquation(f,df,1,0.5e-6)%[xkt]=NewdonToEquation(f,df,1,0.5e-6)%函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6%[xkt]=NewdonToEquation(f,df,1)ifnargin==3eps=0.5e-6;endtic;k=0;while1x=x0-f(x0)./df(x0);k=k+1;ifabs(x-x0)eps||k==30break;endx0=x;endt=toc;ifk==30disp('迭代次数太多。');x='';end2、Newdon迭代法求解非线性方程组functiony=NewdonF(x)%牛顿迭代法解非线性方程组的测试函数%定义是必须定义为列向量y(1,1)=x(1).^2-10*x(1)+x(2).^2+8;y(2,1)=x(1).*x(2).^2+x(1)-10*x(2)+8;return;functiony=NewdonDF(x)%牛顿迭代法解非线性方程组的测试函数的导数y(1,1)=2*x(1)-10;y(1,2)=2*x(2);y(2,1)=x(2).^+1;y(2,2)=2*x(1).*x(2)-10;return;以上两个函数仅供下面程序的测试function[xkt]=NewdonToEquations(f,df,x0,eps)%牛顿迭代法解非线性方程组%[xkt]=NewdonToEquations(f,df,x0,eps)%x:近似解%k:迭代次数%t:运算时间%f:方程组（事先定义）%df:方程组的导数（事先定义）%x0:初始值%eps:误差限%%说明：由于虚参f和df的类型都是函数，使用前需要事先在当前目录下采用函数M文件定义%另外在使用此函数求解非线性方程组时，需要在函数名前加符号“@”，如下所示%%应用举例：%x0=[0,0];eps=0.5e-6;%x=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)%[xk]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)%[xkt]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)%函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6%[xkt]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)ifnargin==3eps=0.5e-6;endtic;k=0;while1x=x0-inv(df(x0))*f(x0);%此处可采用其他方法避免求逆k=k+1;ifnorm(x-x0)eps||k==15break;endx0=x;endt=toc;ifk==15disp('迭代次数太多。');x='';end3、Lagrange插值法提供两个程序，采用了不同的方法functionf=InterpLagrange(x,y,x0)%构造Lagrange插值多项式%此函数中借助向量卷积来求Lagrange基函数，运算速度较快%f=InterpLagrange(x,y,x0)%f:插值多项式或者是插值多项式在x0处的值%x:节点%y:函数值%x0:某一测试点%%调用格式:%f=InterpLagrange(x,y)返回插值多项式%f=InterpLagrange(x,y,x0)返回插值多项式在点x0处的值%举例:%x=[0.320.340.36];y=[0.3145670.3334870.352274];x0=0.33;%f=InterpLagrange(x,y)%f=InterpLagrange(x,y,x0)iflength(x)==length(y)n=length(x);elsedisp('节点个数和函数值个数不同！')f='';return;endp=0;fori=1:nl=y(i);forj=1:nifj==icontinue;end%利用卷积计算Lagrange基函数l=conv(l,[1-x(j)]./(x(i)-x(j)));end%p是一向量，表示插值多项式的系数p=p+l;endifnargin==3f=polyval(p,x0);%计算插值多项式在x0处的值elsef=poly2str(p,'x');%把插值多项式的向量形式转化为插值多项式的符号形式endfunctionf=InterpLagrange2(x,y,x0)%构造Lagrange插值多项式%此函数中借助符号运算来求Lagrange基函数，运算速度较慢，不推荐此种方法%f=InterpLagrange2(x,y,x0)%f:插值多项式或者是插值多项式在x0处的值%x:节点%y:函数值%x0:某一测试点%%调用格式:%f=InterpLagrange2(x,y)返回插值多项式%f=InterpLagrange2(x,y,x0)返回插值多项式在点x0处的值%举例:%x=[0.320.340.36];y=[0.3145670.3334870.352274];x0=0.33;%f=InterpLagrange2(x,y)%f=InterpLagrange2(x,y,x0)iflength(x)==length(y)n=length(x);elsedisp('节点个数和函数值个数不同！')f='';return;endsymst;f=0;fori=1:nl=y(i);forj=1:nifj==icontinue;endl=l*(t-x(j))/(x(i)-x(j));%借助符号运算，计算Lagrange基函数endf=f+l;simplify(f);%化简多项式ifi==nifnargin==3f=subs(f,'t',x0);%计算插值多项式f在点x0处的值elsef=collect(f);%计算插值多项式，展开并合并同类项f=vpa(f,6);%设置多项式系数的有效数字endendend4、Newdon插值法functionf=InterpNewdon(x,y,x0)%Newdon插值多项式%f=InterpNewdon(x,y,x0)%f:插值多项式或者是插值多项式在x0处的值%x:节点%y:函数值%x0:某一测试点%%调用格式%f=InterpNewdon(x,y)返回插值多项式%f=InterpNewdon(x,y,x0)返回插值多项式在x0点的值%应用举例:%x=[12345];y=[14786];x0=6;%f=InterpNewdon(x,y)%f=InterpNewdon(x,y,x0)iflength(x)==length(y)n=length(x);elsedisp('节点个数和函数值个数不同！')f='';return;endA=zeros(n);%初始化差商矩阵fori=1:nA(i,1)=y(i);%差商矩阵的第一列是函数值end%计算差商矩阵%差商矩阵中对角线上的元素为Newdon插值多项式的系数forj=2:nfori=j:nA(i,j)=(A(i,j-1)-A(i-1,j-1))/(x(i)-x(i-j+1));endend%求Newdon插值多项式p=zeros(1,n);fori=1:np1=A(i,i);%差商矩阵对角线上的元素就是Newdon插值多项式的系数forj=1:i-1p1=conv(p1,[1-x(j)]);%计算Newdon插值多项式的基项endp1=[zeros(1,n-i),p1];%向量相加，维数必须相同。把向量的元素补齐p=p+p1;endifnargin==3f=polyval(p,x0);%计算插值多项式在x0处的值elsef=poly2str(p,'x');%把插值多项式的向量形式转化为插值多项式的符号形式end5、基本Guass消去法求解线性方程组functionx=EqtsBasicGuass(A,b)%基本Guass消去法求解线性方程组Ax=b%x=EqtsBasicGuass(A,b)%x:解向量，列向量%A:线性方程组的矩阵%b:列向量%%应用举例:%A=[223;477;-245];b=[3;1;-7];%x=EqtsBasicGuass(A,b)%检查输入参数ifsize(A,1)~=size(b,1)disp('输入参数有误！');x='';return;end%(A|b)A=[Ab];%消去过程n=size(A,1);l=zeros(n);fork=1:n-1fori=k+1:nl(i,k)=A(i,k)/A(k,k);endfori=k+1:nforj=k+1:n+1A(i,j)=A(i,j)-l(i,k)*A(k,j);endforj=1:kA(i,j)=0;endendend%回代过程x=zeros(n,1);x(n)=A(n,n+1)/A(n,n);fori=n-1:-1:1y=0;forj=i+1:ny=y+A(i,j)*x(j);endx(i)=(A(i,n+1)-y)/A(i,i);endreturn;6、三角分解法求解线性方程组functionx=EqtsDoolittleLU(A,b)%Doolittle分解法求解线性方程组Ax=b%x=EqtsDoolittleLU(A,b)%x:解向量，列向量%A:线性方程组的矩阵%b:列向量%%应用举例:%A=[621-1;2410;114-1;-10-13];b=[6;-1;5;-5];%x=EqtsDoolittleLU(A,b)%检查输入参数ifsize(A,1)~=size(b,1)disp('输入参数有误！');x='';return;end%分解%把L和U的元素存储在A的相应位置上n=length(b);fork=1:nforj=k:nz=0;forr=1:k-1z=z+A(k,r)*A(r,j);endA(k,j)=A(k,j)-z;endfori=k+1:nz=0;forr=1:k-1z=z+A(i,r)*A(r,k);endA(i,k)=(A(i,k)-z)/A(k,k);endend%求解x=zeros(size(b));fori=1:nz=0;fork=1:i-1z=z+A(i,k)*x(k);endx(i)=b(i)-z;endfori=n:-1:1z=0;fork=i+1:nz=z+A(i,k)*x(k);endx(i)=(x(i)-z)/A(i,i);endreturn7、追赶法求解三对角线性方程组functionx=EqtsForwardAndBackward(L,D,U,b)%追赶法求解三对角线性方程组Ax=b%x=EqtsForwardAndBackward(L,D,U,b)%x:三对角线性方程组的解%L:三对角矩阵的下对角线，行向量%D:三对角矩阵的对角