目录
6.1数值微分与数值积分
数值微分
- 微积分中,任意函数f(x)在点的导数是通过极限定义的,因此可以利用差商近似代替在点处的导数.
- 数值微分的实现
dx = diff(x) 计算向量x的向前差分,dx(i)=x(i+1)-x(i)
dx = diff(x, n) 计算向量x的n阶向前差分
dx = diff(A, n, dim) 计算矩阵的n阶差分,dim=1时(默认状态),按列计算差分,dim=2,按行计算差分
数值积分
- 常使用牛顿-莱布尼兹公式求函数的积分,或者用数值解法求解定积分的近似值,基本思想是将区间分解为n个子区间,进行求和
- 数值积分的实现
[l, n] = quad(filename, a, b, tol, trace ) 基于自适应辛普森方法
[l, n] = quadl(filename, a, b, tol, trace ) 基于自适应辛普森方法
其中,filename是被积函数名; a和b是定积分的上限和下限(积分限必须是有限的),tol用来控制积分情况,默认时取,trace控制是否展现积分过程,若取非0则展现积分过程,取0则不展现(默认),返回参数I即为定积分的值,n为被积函数的调用次数.
l = integral(filename, a, b ) 基于全局自适应积分方法
其中,,filename是被积函数名;a和b是定积分的上限和下限(积分限可以为无穷大),返回参数I即为定积分的值
[l, err] = quadgk(filename, a, b) 基于自适应高斯-克朗罗德方法
其中,err返回近似误差范围,其他参数的含义和用法与quad函数相同.积分上下限可以实无穷大(-Inf或Inf),也可以是复数.如果积分上下限是复数,则quadgk函数在复平面上求积分.
l = trapz(x, y) 基于梯形积分法
其中,向量x,y定义函数关系y=f(x)
- 多重定积分的数值求解
l = integral2(filename, a, b, c, d)
l = quad2d(filename, a, b, c, d)
l = dblquad(filename, a, b, c, d, tol)
l = integral3(filename, a, b, c, d, e, f)
l = triplequad(filename, a, b, c, d, e, f, tol)
6.2线性方程求解
直接法
- 利用左除运算符的直接解法
Ax=b,因此x=A\b
若右端项为n*m的矩阵,则x=A\b可同时获得稀疏矩阵A相同的m个线性方程组的数值解x,x为n*m的矩阵.即x(:,j)=A\b(:,j),j=1,2,...,m
A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];
b=[13,-9,6,0]';
x=A\b
- 利用矩阵分解求解线性方程组
主要分为以下三种分解方式
- LU分解(以此为例)
- QR分解
- Cholesky分解
基本思想如下:
分解函数
[L,U]=lu(A) 产生一个上三角阵U和一个变换形式的下三角阵L,使之满足A=LU,注意,这里的矩阵A必须是方阵,最后x=U\(L\b)
[L,U,P]=lu(A) 产生一个上三角阵U和一个下三角阵以及一个置换矩阵P,使之满足PA=LU,最后x=U\(L\P*b)
A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];
b=[13,-9,6,0]';
[L,U]=lu(A);
x=U\(L\b)
迭代法
- 雅可比迭代法
如果x收敛了,则就是所解的值.
function [y,n]=jacobi(A,b,x0,ep)
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
B=D\(L+U);
f=D\b;
y=B*x0+f;
n=1;
while norm(y-x0)>=ep
x0=y;
y=B*x0+f;
n=n+1;
end
x0为设置的初值,n为迭代次数,下同
- 高斯-赛德尔迭代法(精度更高一些)
function [y,n]=gauseidel(A,b,x0,ep)
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
B=(D-L)\U;
f=(D-L)\b;
y=B*x0+f;
n=1;
while norm(y-x0)>=ep
x0=y;
y=B*x0+f;
n=n+1;
end
6.3非线性方程求解与函数极值计算
非线性方程求解
x=fzero(filename,x0) 其中filename是待求解方程左端的函数表达式,x0是初始值
x=fsolve(filename,x0,option) option为优化参数,可设置为optimset('Display','off'),可以求解非线性方程组
f=@(x) [sin(x(1))+x(2)+x(3)^2*exp(x(1)),x(1)+x(2)+x(3),x(1)*x(2)*x(3)];
f([1,1,1])
x=fsolve(f,[1,1,1],optimset('Display','off'))
f(x)
函数极值计算
- 无约束最优化问题
[xmin, fmin] = fminbnd(filename, x1, x2, option)
[xmin, fmin] = fminsearch(filename, x0, option)
[xmin, fmin] = fminunc(filename, x0, option)
其中x1,x2为研究区间的左右边界,x0是表示极值点的初值
- 有约束最优化问题
- 线性不等式约束
- 线性等式约束
- 非线性不等式约束
- 非线性等式约束
- x的上界和下界
[xmin, fmin] = fmincon(filename, x0, A, b, Aep, bep, Lbnd, Ubnd, NonF, option)
其中xmin, fmin, fiename, x0, option的含义与求最小值函数相同.其余参数为约束条件,包括线性不等式约束,线性等式约束,x的下界和上界以及定义非线性约束的函数.如果某个约束不存在,则用空矩阵表示.
6.4常微分方程数值求解
一般概念
- 寻找函数y满足y'=f(t, y)
求解函数
[t, y] = solver(filename, tspan, y0, option)
- 其中,t和y分别给出时间向量和相应的数值解.solver为求常微分返程数值解的函数.filename是定义f(t,y)的函数名,该函数必须返回一个列向量.tspan形式为[t0, tf],表示求解区间,y0是初始状态列向量,option是可选参数,用于设置求解属性.
- solver函数有统一的命名方式,odennxx.nn表示左右方法的阶数,常用ode45
f=@(t,y) (y^2-t-2)/4/(t+1);
[t,y]=ode23(f,[0,10],2);
y1=sqrt(t+1)+1;
plot(t,y,'b:',t,y1,'r');
刚性问题
- 解的分量有的变化很快,有的变化很慢,而且相差悬殊,这就是所谓的刚性问题(Stiff)
求解函数ode15s,使用方法和上面的一样