1、Newdon迭代法求解非线性方程function [x k t]=NewdonToEquation(f,df,x0,eps)%牛顿迭代法解线性方程%[x k t]=NewdonToEquation(f,df,x0,eps)%x:近似解%k:迭代次数%t:运算时间%f:原函数，定义为内联函数�:函数的倒数，定义为内联函数%x0:初始值%eps:误差限%%应用举例：%f=inline('x^3+4*x^2-10');�=inline('3*x^2+8*x');%x=NewdonToEquation(f,df,1,0.5e-6)%[x k]=NewdonToEquation(f,df,1,0.5e-6)%[x k t]=NewdonToEquation(f,df,1,0.5e-6)%函数的最后一个参数也可以不写。

默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquation(f,df,1)if nargin==3eps="0".5e-6;endtic;k=0;while 1x="x0-f"(x0)./df(x0);k="k"+1;if abs(x-x0) < eps || k >30break;endx0=x;endt=toc;if k >= 30disp('迭代次数太多。

');x="0";t="0";end2、Newdon迭代法求解非线性方程组function y="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;function y="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 [x k t]=NewdonToEquations(f,df,x0,eps)%牛顿迭代法解非线性方程组%[x k t]=NewdonToEquations(f,df,x0,eps)%x:近似解%k:迭代次数%t:运算时间%f:方程组（事先定义）�:方程组的导数（事先定义）%x0:初始值%eps:误差限%%说明：由于虚参f和df的类型都是函数，使用前需要事先在当前目录下采用函数M文件定义% 另外在使用此函数求解非线性方程组时，需要在函数名前加符号“@”，如下所示%%应用举例：%x0=[0,0];eps=0.5e-6;%x=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)%[x k]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)%[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)%函数的最后一个参数也可以不写。

默认情况下，eps=0.5e-6%[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)if nargin==3eps="0".5e-6;endtic;k=0;while 1x="x0-inv"(df(x0))*f(x0);%此处可采用其他方法避免求逆k="k"+1;if norm(x-x0) < eps || k > 15break;endx0=x;endt=toc;if k >= 15disp('迭代次数太多。

');x="zeros"(size(x0));t="0";end3、Lagrange插值法提供两个程序，采用了不同的方法function f="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.32 0.34 0.36];y=[0.314567 0.333487 0.352274];x0=0.33;%f=InterpLagrange(x,y)%f=InterpLagrange(x,y,x0)if length(x)==length(y)n="length"(x);elsedisp('节点个数和函数值个数不同！')f=' ';return;endp=0;for i="1:n"l="y"(i);for j="1:n"if j==icontinue;end%利用卷积计算Lagrange基函数l=conv(l,[1 -x(j)]./(x(i)-x(j)));end%p是一向量，表示插值多项式的系数p="p"+l;endif nargin==3f="polyval"(p,x0);%计算插值多项式在x0处的值elsef="poly2str"(p,'x');%把插值多项式的向量形式转化为插值多项式的符号形式endfunction f="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.32 0.34 0.36];y=[0.314567 0.333487 0.352274];x0=0.33; %f=InterpLagrange2(x,y)%f=InterpLagrange2(x,y,x0)if length(x)==length(y)n="length"(x);elsedisp('节点个数和函数值个数不同！')f=' ';return;endsyms t;f=0;for i="1:n"l="y"(i);for j="1:n"if j==icontinue;endl=l*(t-x(j))/(x(i)-x(j));%借助符号运算，计算Lagrange基函数endf="f"+l;simplify(f);%化简多项式if i==nif nargin==3f="subs"(f,'t',x0);%计算插值多项式f在点x0处的值elsef="collect"(f);%计算插值多项式，展开并合并同类项f="vpa"(f,6);%设置多项式系数的有效数字endendend4、Newdon插值法function f="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=[1 2 3 4 5];y=[1 4 7 8 6];x0=6;%f=InterpNewdon(x,y)%f=InterpNewdon(x,y,x0)if length(x)==length(y)n="length"(x);elsedisp('节点个数和函数值个数不同！')f=' ';return;endA=zeros(n);%初始化差商矩阵for i="1:n"A(i,1)=y(i);%差商矩阵的第一列是函数值end%计算差商矩阵%差商矩阵中对角线上的元素为Newdon插值多项式的系数for j="2:n"for i="j:n"A(i,j)=(A(i,j-1)-A(i-1,j-1))/(x(i)-x(i-j+1));endend%求Newdon插值多项式p=zeros(1,n);for i="1:n"p1=A(i,i);%差商矩阵对角线上的元素就是Newdon插值多项式的系数 for j="1:i-1"p1=conv(p1,[1 -x(j)]);%计算Newdon插值多项式的基项endp1=[zeros(1,n-i),p1];%向量相加，维数必须相同。

把向量的元素补齐p="p"+p1;endif nargin==3f="polyval"(p,x0);%计算插值多项式在x0处的值elsef="poly2str"(p,'x');%把插值多项式的向量形式转化为插值多项式的符号形式end5、基本Guass消去法求解线性方程组function x="EqtsBasicGuass"(A,b)%基本Guass消去法求解线性方程组Ax=b%x=EqtsBasicGuass(A,b)%x:解向量，列向量%A:线性方程组的矩阵%b:列向量%%应用举例:%A=[2 2 3;4 7 7;-2 4 5]; b=[3;1;-7];%x=EqtsBasicGuass(A,b)%检查输入参数if size(A,1) ~= size(b,1)disp('输入参数有误！');x=' ';return;end%(A|b)A=[A b];%消去过程n=size(A,1);l=zeros(n);for k="1:n-1"for i="k"+1:nl(i,k)=A(i,k)/A(k,k);endfor i="k"+1:nfor j="k"+1:n+1A(i,j)=A(i,j)-l(i,k)*A(k,j);endfor j="1:k"A(i,j)=0;endendend%回代过程x=zeros(n,1);x(n)=A(n,n+1)/A(n,n);for i="n-1:-1:1"y="0";for j="i"+1:ny=y+A(i,j)*x(j);endx(i)=(A(i,n+1)-y)/A(i,i);endreturn;6、三角分解法求解线性方程组function x="EqtsDoolittleLU"(A,b)%Doolittle分解法求解线性方程组Ax=b%x=EqtsDoolittleLU(A,b)%x:解向量，列向量%A:线性方程组的矩阵%b:列向量%%应用举例:%A=[6 2 1 -1;2 4 1 0;1 1 4 -1;-1 0 -1 3];b=[6;-1;5;-5]; %x=EqtsDoolittleLU(A,b)%检查输入参数if size(A,1) ~= size(b,1)disp('输入参数有误！');x=' ';return;end%分解%把L和U的元素存储在A的相应位置上n=length(b);for k="1:n"for j="k:n"z=0;for r="1:k-1"z="z"+A(k,r)*A(r,j);endA(k,j)=A(k,j)-z;endfor i="k"+1:nz=0;for r="1:k-1"z="z"+A(i,r)*A(r,k);endA(i,k)=(A(i,k)-z)/A(k,k);endend%求解x=zeros(size(b));for i="1:n"z="0";for k="1:i-1"z=z+A(i,k)*x(k);endx(i)=b(i)-z;endfor i="n:-1:1"z="0";for k="i"+1:nz=z+A(i,k)*x(k);endx(i)=(x(i)-z)/A(i,i);endreturn7、追赶法求解三对角线性方程组function x="EqtsForwardAndBackward"(L,D,U,b)%追赶法求解三对角线性方程组Ax=b%x=EqtsForwardAndBackward(L,D,U,b)%x:三对角线性方程组的解%L:三对角矩阵的下对角线，行向量%D:三对角矩阵的对角线，行向量%U:三对角矩阵的上对角线，行向量%b:线性方程组Ax=b中的b，列向量%%应用举例:%L=[-1 -2 -3];D=[2 3 4 5];U=[-1 -2 -3];b=[6 1 -2 1]'; %x=EqtsForwardAndBackward(L,D,U,b)%检查参数的输入是否正确n=length(D);m=length(b);n1=length(L);n2=length(U);if n-n1 ~= 1 || n-n2 ~= 1 || n ~= mdisp('输入参数有误！')x=' ';return;end%追的过程for i="2:n"L(i-1)=L(i-1)/D(i-1);D(i)=D(i)-L(i-1)*U(i-1);endx=zeros(n,1);x(1)=b(1);for i="2:n"x(i)=b(i)-L(i-1)*x(i-1);end%赶的过程x(n)=x(n)/D(n);for i="n-1:-1:1"x(i)=(x(i)-U(i)*x(i+1))/D(i);endreturn;8、主元素的Guass消去法求解线性方程组function x="EqtsClmnPrimElemGuass"(A,b)%主元素的Guass消去法求解线性方程组Ax=b%x=EqtsClmnPrimElemGuass(A,b)%x:解向量，列向量%A:线性方程组的矩阵%b:列向量%%应用举例:%A=[-0.002 2 2;1 0.78125 0;3.996 5.5625 4];%b=[0.4;1.3816;7.4178];%x=EqtsClmnPrimElemGuass(A,b)%检查输入参数if size(A,1) ~= size(b,1)disp('输入参数有误！');x=' ';return;end%(A|b)A=[A b];%消去过程n=size(A,1);l=zeros(n);for k="1:n-1"%换行[a idx1]=max(abs(A(k:n,k)));%寻找绝对值最大的元素的下标[b idx2]=min(abs(A(k:n,k)));%寻找绝对值最小的元素的下标 idx1=idx1+k-1;idx2=idx2+k-1;for j="1:n"+1c=A(idx1,j);A(idx1,j)=A(idx2,j);A(idx2,j)=c;endfor i="k"+1:nl(i,k)=A(i,k)/A(k,k);endfor i="k"+1:nfor j="k"+1:n+1A(i,j)=A(i,j)-l(i,k)*A(k,j);endfor j="1:k"A(i,j)=0;endendend%回代过程x=zeros(n,1);x(n)=A(n,n+1)/A(n,n);for i="n-1:-1:1"y="0";for j="i"+1:ny=y+A(i,j)*x(j);endx(i)=(A(i,n+1)-y)/A(i,i);endreturn;9、Euler法求解常微分方程function outXY="ODEEuler"(f,x0,y0,h,PointNum)%简单欧拉法求解常微分方程dy/dx=f%outXY=ODEEuler(f,x0,y0,h,PointNum)%outXY:所取点的横纵坐标。