Matlab 数值分析计算汇集

%追赶法求解线性方程组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);
    end
end
x(n)=y(n);
for j=n-1:-1:1
    x(j)=y(j)-x(j+1)*u(j);
end

clc
clear
A=[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)

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
    end
end

x=[11,12,13];
y=[2.3979,2.4849,2.5649];
xh=11.75;
yh=lagrange(x,y,xh)

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))';  %求差商表  （注意这里有一个 ’ 符号，与差商表不一样的地方）
end
q=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)

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)

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)

%高斯-赛德尔迭代法：
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;
    end
end

A=[2,-1,0;-1,3,-1;0,-1,2];
b=[1;8;-5];
tol=1e-4;
[x,iter]=gs(A,b,tol)

%牛顿法求非线性方程的根：
%   输入：fun--非线性函数；dfun--非线性函数导数；x0--初始值；tol--精度；
%   输出：x--非线性方程数值根
function [x,iter]=newton(fun,dfun,x0,tol)
format long
iter=1;
x=x0;
while iter<500
  x=x-feval(fun,x)/feval(dfun,x);
  if abs(feval(fun,x))<tol
    break;
  end
  iter=iter+1;
end

function y=fun3(x)
y=x*cos(x)+2;%

function y=dfun3(x)
y=cos(x)-x*sin(x);

x=newton(@fun3,@dfun3,2,1e-3)

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;
end
t=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)

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;
end
t=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)

function y=rk4(fun,a,b,y0,n)
h=(b-a)/n;
y(1)=y0;
for k=1:n
  x=a+(k-1)*h;
  k1=h*feval(fun,x,y(k));
  k2=h*feval(fun,x+h/2,y(k)+k1/2);
  k3=h*feval(fun,x+h/2,y(k)+k2/2);
  k4=h*feval(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)

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)