数值分析中常用的matlab程序

1.

% 最小二乘法拟合数据点方法1：

% 左除右除：xA=B ==> x=B/A | Ax=B ==> x=B\A

% A 表示拟合函数的组合，如：多项式插值，A=[1,x,x.^2,...,x.^n],表示拟合函数为

% 多项式：s(x)=a0+a1*x+...+an*x^n;又如：A=[log(x),cos(x),exp(x)]

%则表示拟合函数为s(x)=a0*ln(x)+a1*cos(x)+a2*exp(x)

% 法方程为：A'*A*z=A'*y ==> A*z=y ==> z=A\y z=(a0,...,an)'

% Date: 2012-1-1

clear;

clc;

x=[0 0.25 0.5 0.75 1]';

y=[1 1.284 1.6487 2.1170 2.7183]';%索要拟合的数据点

x1=ones(siz e(x),1);

A=[x1 x x.^2];%拟合函数

Z=A\y %A中每列函数的参数

% 最小二乘法拟合数据方法2：采用polyfit函数

% Date:2012-1-1

clear;

clc;

x=[0 0.25 0.5 0.75 1]';

y=[1 1.284 1.6487 2.1170 2.7183]';%索要拟合的数据点

p=polyfit(x,y,2) %polyf it(x,y,n),(x,y)为数据点坐标，n为拟合多项式阶数，p 为

% 所求拟合多项式的幂次从高到低排列的系数。

2.

% 复合梯形公式计算积分值

% 输入：fun--积分函数；a，b--积分区间；n--区间等分数

% 输出：I--数值积分结果

% 调用格式(ex)：re=ftrapz(@fun1,0,1,10)

% 2012-1-1

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

h=(b-a)/n;%区间等分

x=linspace(a,b,n+1);%将a到b的区间等分成(n+1)-1个区间，数据点有(n+!)个

y=feval(fun,x);

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

%积分原函数

% Date：2012-1-1

function y=fun1(x)

y=exp(-x);

3.

% 复合simpson公式求积分

% 输入：fun--积分函数；a，b--积分区间；n--区间等分数% 输出：I--数值积分结果

% 调用格式(ex)：re=fsimpson(@fun1,0,1,10)

% 2012-1-1

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));

4.

% 两点GS-Legendre公式求积分

% 输入：fun--积分函数；a，b--积分区间；

% 输出：I--数值积分结果

% 调用格式(ex)：re=GSLege(@fun2,0,1)

% 2012-1-1

function I=GSLege(fun,a,b)

%将区间[a.b]通过变量替换x=(a+b)/2-(b-a)/2*t变到[-1,1]

%其中t取GS-Lege的高斯点

m1=f eval(fun,(a+b)/2+(b-a)/2*(-1/sqrt(3)));

m2=f eval(fun,(a+b)/2+(b-a)/2*(1/sqrt(3)));

I=(b-a)/2*(m1+m2);

% GS-Legendre公式的积分函数

% Date:2012-1-1

function y=fun2(x)

y=sin(x)/x;

5.

% 追赶法求解线性方程组Ax=b，其中A是三对角方阵

%function x=tridiagsolver(A,b)

clear;

clc;

A=[2 -1 0 0;-1 3 -2 0;0 -2 4 -3;0 0 -3 5];%三对角矩阵，线性方程组系数矩阵b=[6,1,-2,1]';%

[n,n]=size(A);

fori=1:n

if(i<2)

l(i)=A(i,i);

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

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

elseif i

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);

u(i)=A(i,i+1)/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

end

x(n)=y(n);

for j=n-1:-1:1

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

end

x

6.

% SOR迭代求解非线性方程组Ax=b

% 输入：A--系数矩阵；b--；omega--松弛因子（0~2);tol--精度

% 输出：x--方程的解向量;iter--迭代次数;

% 调用格式(ex)：[x,iter=sor(A,b,1.1,1e-4)

% 2012-1-1

%function [x,iter]=sor(A,b,omega,tol)

%{%}

clear;clc;

A=[2 -1 0;-1 3 -1;0 -1 2];

b=[1 8 -5]';

omega=1.1;

tol=1e-4;

D=diag(diag(A));%diag(A)返回的是A的对角元组成的列向量；diag(b)返回的是以列向量b为对角元的方阵；

L=D-tril(A);%tril(A)返回的是矩阵A的下三角矩阵；

U=D-triu(A);%triu(A)返回的是矩阵A的上三角矩阵；

x=zeros(size(b));%给定迭代初始值零向量

iter=1;

whileiter<500

x=(D-omega*L)\(omega*b+(1-omega)*D*x+omega*U*x);%迭代格式

error=norm(b-A*x)/norm(b);%收敛条件

if error

break;

end

iter=iter+1;

end

ifiter>= 500

fprintf('root not found!');

end

7.

% 牛顿法求解非线性方程的根

% 输入：fun--非线性函数；dfun--非线性函数导数；x0--初始值；tol--精度；

% 输出：x--非线性方程数值根

% 调用格式(ex)：x=newton(@fun3,@dfun3,6,1e-3)

% 2012-1-1

function x=newton(fun,dfun,x0,tol)

iter=1;

if abs(f eval(fun,x0))

x=x0; return;

end

whileiter<500

ifiter==1

x=x0-feval(fun,x0)/f eval(dfun,x0);

if abs(f eval(fun,x))

break;

end

else

x=x-feval(fun,x)/feval(dfun,x);

if abs(f eval(fun,x))

break;

end

end

iter=iter+1;

end

ifiter==500

fprintf('not s uccessf ul!');

x=NaN;

end

% newton的函数文件

% Date：2012-1-1

function y=fun3(x)

y=3.*x.^3-8.*x.^2-8.*x-11;

% newton的导函数文件

% Date：2012-1-1

function y=dfun3(x)

y=9.*x.^2-16.*x-8;

8.

% 两点割线法求解非线性方程的根

% 输入：fun--非线性函数；a,b--两个初始值；tol--精度；

% 输出：x--非线性方程数值根

% 调用格式(ex)：x=gexian(@fun3,3,6,1e-3)

% 2012-1-1

function x=gexian(fun,a,b,tol)

iter=1;

xk1=a;

xk2=b;

whileiter<500

xk3=xk2-feval(fun,xk2)*(xk2-xk1)/(feval(fun,xk2)-f eval(fun,xk1)); if abs(f eval(fun,xk3))

break;

else

xk1=xk2;

xk2=xk3;

end

iter=iter+1;

end

ifiter==500

fprintf('not s uccessf ul!');

x=NaN;

else

x=xk3;

end

9.

% 乘幂法求矩阵的按模最大特征值及其特征向量

% 输入：A--要求的矩阵；v0--初始非零向量；tol--精度；% 输出：x--特征向量；lam--按模最大特征值

% 调用格式(ex)：[lam,x]=power(A,v0,1e-2)

% 2012-1-1

function [lam,x]=matrixpower(A,v0,tol)

v1=A*v0;

v2=A*v1;

sum=0;

p=0;

fori=1:siz e(v1)

if v1(i)~=0

sum=sum+v2(i)/v1(i);

p=p+1;

end

end

lam0=sum/p;

iter=2;

vk1=v2;

whileiter<500

vk2=A*vk1;

sum=0;

p=0;

fori=1:siz e(vk1)

if vk1(i)~=0

sum=sum+vk2(i)/vk1(i);

p=p+1;

end

end

lam=sum/p;

if abs(lam-lam0)

break;

else

lam0=lam;

vk1=vk2;

end

end

ifiter==500

fprintf('not s uccessf ul!'); lam=NaN;

else

x=vk2;

end

9_2

% 改进乘幂法求矩阵的按模最大特征值及其特征向量

% 输入：A--要求的矩阵；v0--初始非零向量；tol--精度；

% 输出：x--特征向量；lam--按模最大特征值

% 调用格式(ex)：[lam,x]=pMatrixPower(A,v0,1e-2)

% Date:2012-1-2

function [lam,x]=pMatrixPower(A,v0,tol)

[ty,ti]=max(abs(v0));%返回v0中元素绝对值最大的元素值与下标,ti为下标lam0=v0(ti);

u0=v0/lam0;

iter=1;

whileiter<500

v1=A*u0;

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

lam1=v1(ti);

u0=v1/lam1;

if abs(lam0-lam1)

break;

else

lam0=lam1;

iter=iter+1;

end

end

ifiter>=500

fprintf('not s uccessf ul!');

lam=NaN;

else

lam=lam1;

x=u0;

end

10.

% 反幂法求矩阵的按模最小特征值及其特征向量

% 输入：A--要求的矩阵；v0--初始非零向量；tol--精度；

% 输出：x--特征向量；lam--按模最大特征值

% 调用格式(ex)：[lam,x]=invM atrixPower(A,v0,1e-2)

% Date:2012-1-2

function [lam,x]=invM atrixPower(A,v0,tol)

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

lam0=v0(ti);

u0=v0/lam0;

iter=1;

whileiter<500

v1=A\u0;

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

lam1=v1(ti);

u0=v1/lam1;

if abs(1/lam0-1/lam1)

break;

else

lam0=lam1;

iter=iter+1;

end

end

ifiter>=500

fprintf('not s uccessf ul!');

lam=NaN;

else

lam=1/lam1;

x=u0;

end

11.

% 欧拉方法求解一阶常微分方程初值问题

% 输入：fun--一阶常微分函数；a,b--求解区间；y0--函数在a点值y(a)；n--所分区间数;

% 输出：y--常微分方程在区间[a,b]上各点的数值解；

% 调用格式(ex)：y=odeEuler(@fun4,0,1,1,10)

% Date:2012-1-2

function y=odeEuler(fun,a,b,y0,n)

h=(b-a)/n;

y(1)=y0;

x=a:h:b;

fori=1:n

y(i+1)=y(i)+h*feval(fun,x(i),y(i));

end

% 常微分方程fun4=0

% Date：2012-1-2

function re=fun4(x,y)

re=y-2*x/y;

12.

% 改进欧拉公式（预估--校正）方法求解一阶常微分方程初值问题

% 输入：fun--一阶常微分函数；a,b--求解区间；y0--函数在a点值y(a)；n--所分区间数;

% 输出：y--常微分方程在区间[a,b]上各点的数值解；

% 调用格式(ex)：y=pOdeEuler(@fun4,0,1,1,10)

% Date:2012-1-2

function y=pOdeEuler(fun,a,b,y0,n)

h=(b-a)/n;

y(1)=y0;

x=a:h:b; fori=1:n

yp=y(i)+h*feval(fun,x(i),y(i));

yc=y(i)+h*feval(fun,x(i+1),yp);

y(i+1)=0.5*(yp+yc);

end

13.

% 梯形方法求解一阶常微分方程初值问题

% 输入：fun--一阶常微分函数；a,b--求解区间；y0--函数在a点值y(a)；n--所分区间数;

% 输出：y--常微分方程在区间[a,b]上各点的数值解；

% 调用格式(ex)：y=trapezium(@fun4,0,1,1,10)

% Date:2012-1-2

function y=trapez ium(fun,a,b,y0,n)

h=(b-a)/n;

y(1)=y0;

x=a:h:b;

tol=1e-6;

fori=1:n

%用不动点迭代的方法求解非线性方程：

%y(i+1)=y(i)+h*feval(fun,x(i),y(i))/2+h*feval(fun,x(i),y(i+1))/2;

iter=1;

yy0=1+(i-1)*h;%迭代初始值

whileiter<500

yy1=y(i)+h*feval(fun,x(i),y(i))/2+h*f eval(fun,x(i)+h,yy0)/2;

if abs(yy1-yy0)

break;

else

yy0=yy1;

iter=iter+1;

end

end

ifiter>=500

printf('not s uccessful!');

y=NaN;

return;

else

y(i+1)=yy1;

end

end

14.

% 标准四阶四段龙格-库塔方法求解一阶常微分方程初值问题

% 输入：fun--一阶常微分函数；a,b--求解区间；y0--函数在a点值y(a)；n--所分区间数;

% 输出：y--常微分方程在区间[a,b]上各点的数值解；

% 调用格式(ex)：y=longekuta(@fun4,0,1,1,10)

% Date:2012-1-2

function y=longekuta(fun,a,b,y0,n)

h=(b-a)/n;

y(1)=y0;

for k=2:n+1

x=a+(k-2)*h;

k1=h*feval(fun,x,y(k-1));

k2=h*feval(fun,x+h/2,y(k-1)+k1/2);

k3=h*feval(fun,x+h/2,y(k-1)+k2/2);

k4=h*feval(fun,x+h,y(k-1)+k3);

y(k)=y(k-1)+(k1+2*k2+2*k3+k4)/6;

end

15.插值

% 拉格朗日插值

% 输入：x,y--插值数据点（x，y均为行向量）；xh--要插值的点; % 输出：yh--插值结果；

% 调用格式(ex)：yh=lagrange(x,y,xh)

% Date:2012-1-2

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);

fori=1:n

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

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

%prod对输入的一个向量返回其所有分量的乘积

end

%拉格朗日调用

clear;

clc;

x=[11 12];

y=[2.3979 2.4849];

xh=11.75;

yh=lagrange(x,y,xh)

15_2

% 牛顿插值

% 输入：x,y--插值数据点（x，y均为行向量）；xh--要插值的点; % 输出：yh--插值结果；

% 调用格式(ex)：yh=newtonPol(x,y,xh)

% Date:2012-1-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))';%求差商表end

q=p(1,2:n+1)';%求牛顿法的系数--取第一行

yh=0;

m=1;

yh=q(1);

fori=2:n

m=q(i);

for j=2:i

m=m*(xh-x(j-1));%求牛顿法中各多项式值(xh-x0)…(xh-x n-1) end

yh=yh+m;%求和

end

%牛顿插值调用

clear;

clc;

x=[11 12 13];

y=[2.3979 2.4849 2.5649];

xh=11.75;

yh=newtonPol(x,y,xh)