首页 > 其他分享 >Matlab 数值分析计算汇集

Matlab 数值分析计算汇集

时间:2022-10-12 10:09:08浏览次数:55  
标签:function yh xh 数值 汇集 lam1 Matlab ti fun

分享一下数值分析经常遇到的算法,代码有点多;算法原理之类的网上均可以找到,本文只给出对应的代码实现。

Matlab 数值分析计算汇集_迭代


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)



标签:function,yh,xh,数值,汇集,lam1,Matlab,ti,fun
From: https://blog.51cto.com/domi/5748876

相关文章