数值分析重要算法的matlab程序

数值分析重要算法的matlab程序

插值多项式:拉格朗日插值

function yh=lagrange(x,y,xh)

n=length(x);

m=length(xh);

yh=zeros(1,m);

c1=ones(n-1,1);

c2=ones(1,m);

for i=1:n

xp=x([1:i-1 i+1:n]);

yh=yh+y(i)*prod((c1*xh-xp'*c2)./(x(i)-xp'*c2));

end

插值多项式:牛顿插值（以数值实验3.2）为例

function y=ex32

n = 21;

x = linspace(-5,5,n)';

h = (5-(-5))/(n-1);

y = 1./(1+x.^2);

% form the differences table

for j = 2:n,

y(1:n+1-j,j) = diff(y(1:n+2-j,j-1))./(x(j:n)-x(1:n+1-j));

end

% newton coeff

y = y(1,:);

pz = [ ];

v = linspace(-5,5,80);

for t = v,

z = y(n);

for j = n-1:-1:1,

z = z * (t-x(j))+y(j);

end

pz=[pz z];

end

plot(v,pz,'r+-',v,1./(1+v.^2),'g--');

数值积分：梯形求积公式求积分

function I=ftrapz(fun,a,b,n)

h=(b-a)/n;

x=linspace(a,b,n+1);

y=feval(fun,x);

I=h*(0.5*y(1)+sum(y(2:n))+0.5*y(n+1));

数值积分：抛物型求积公式求积分

function I=fsimpsion(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 x=tridiagsolver(a,b)

[n,n]=size(a);

l=zeros(1,n);

y=zeros(1,n);

u=zeros(1,n);

for i=1:n

if (i==1)

l(i)=a(i,i);

y(i)=b(i)/l(i);

else

l(i)=a(i,i)-a(i,i-1)*u(i-1);

y(i)=(b(i)-y(i-1)*a(i,i-1))/l(i);

end

if (i

u(i)=a(i,i+1)/l(i);

end

end

x(n)=y(n);

for j=n-1:-1:1

x(j)=y(j)-u(j)*x(j+1);

end

高斯-赛德尔迭代解线性方程组

function [x,iter]=gs(A,b,tol)

D=diag(diag(A));

L=D-tril(A);

U=D-triu(A);

x=zeros(size(b));

for iter=1:1000;

x=(D-L)\(b+U*x);

error=norm(b-A*x)/norm(b);

if(error

break;

end

end

牛顿法求方程的根

function [x, it, convg]=newton(x0, f, g, maxit, tol)

%find the zero of function f, with gradient g provided %Usage: [x, it, convg] = newton(x0, f, g, maxit, tol) if nargin<5, tol = 1e-5;

if nargin<4, maxit=100;

end

end

x=x0;

fx=feval(f,x);

convg = 0;

it=1;

while ~convg,

it=it+1;

if norm(fx)

fprintf('Newton Iteration successes!!\n');

convg=1;

return;

end

d=-feval(g,x)\fx;

lambda=1;

lsdone=0;

while ~lsdone,

xn=x+lambda*d';

fn=feval(f,xn);

if abs(fn)

lsdone=1;

else

lambda=1/2*lambda;

if lambda<=eps,

convg=-1;

error('line search fails!!');

end

end

end

x=xn;

fx=fn;

if it > maxit,

convg=0;

error('Newton method needs more iterations!!');

end

end

龙格-库塔方法求解微分方程数值解

function [x,y]=rk4(ef,tspan,y0,n)

y=zeros(1,n+1);

y(1)=y0;

a=tspan(1);

b=tspan(2);

h=(b-a)/n;

x=a:h:b;

c1=[1 2 2 1]'/6;

for i=1:n,

k(1)=h*feval(ef,x(i),y(i));

k(2)=h*feval(ef,x(i)+0.5*h,y(i)+0.5*k(1));

k(3)=h*feval(ef,x(i)+0.5*h,y(i)+0.5*k(2));

k(4)=h*feval(ef,x(i)+h,y(i)+k(3));

y(i+1)=y(:,i)+k*c1;

end

乘幂法求特征值和特征向量

function [t,y]=eigIPower(a,xinit,ep)

v0=xinit;

[tv,ti]=max(abs(v0));

lam0=v0(ti);

u0=v0/lam0;