✅作者简介:热爱科研的算法开发者,Python、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
调用程序
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)
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
end
end
调用程序
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))'; %求差商表 (注意这里有一个 ’ 符号,与差商表不一样的地方)
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)
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
主函数(调用程序)
A=[2,-1,0;-1,3,-1;0,-1,2];
b=[1;8;-5];
tol=1e-4;
[x,iter]=gs(A,b,tol)
5、非线性方程求根
%牛顿法求非线性方程的根:
% 输入:fu服务器托管网n--非线性函数;dfun--非线性函数导数;x0--初始值;tol--精度;
% 输出:x--非线性方程数值根
function [x,iter]=newton(fun,dfun,x0,tol)
format long
iter=1;
x=x0;
while iter
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;
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)
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=Au0;
[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)
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;
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)
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)
服务器托管,北京服务器托管,服务器租用 http://www.fwqtg.net
机房租用,北京机房租用,IDC机房托管, http://www.fwqtg.net