6. 数值微积分与方程求解

A 数值微分与数值积分

A.a数值微分(diff)

<1>数值差分与差商
任意函数 f ( x ) f(x) f(x) x 0 x0 x0点的导数是通过极限定义的:
在这里插入图片描述
在这里插入图片描述

在这里插入图片描述

<2>数值微分的实现
MATLAB提供了求向前差分的函数diff,其调用格式有三种:

  • dx=diff(x):计算向量x的一阶向前差分,dx(i)=x(i+1)-x(i),i=1,2,...,n-1
  • dx=diff(x, n) :计算向量x的n阶向前差分。例如:diff(x, 2)=diff(diff(x))
  • dx=diff(A,n,dim):计算矩阵A的n阶差分,dim=1时(默认状态),按列计算差分;dim=2时,按行计算差分。

注意:diff函数计算的是向量元素间的差分,故差分向量元素的个数比原向量少了一个。同样,对于矩阵来说,差分后的矩阵比原矩阵少了一行或一列。

计算差分之后,可以计算 f ( x ) f(x) fx在某点的差商,计算 f ′ ( x ) f'(x) f(x)的近似值。
例子:

%1*5000的矩阵
x = [0,sort(2*pi*rand(1,5000)),2*pi];
y = sin(x);
f3 = diff(y)./diff(x);
f2 = cos(x(1:end - 1));
plot(x(1:end-1),f3,x(1:end-1),f2);
d = norm(f3-f2)

在这里插入图片描述

A.b 数值积分

<1>数值积分基本原理
牛顿-莱布尼茨公式:
在这里插入图片描述
按积分区间[a, b]分成n个子区间 [ x i , x i + 1 ] , i = 1 , 2 , . . . , n [x_i, x_{i+1}],i=1, 2,...,n [xi,xi+1]i=1,2,...,n,其中 x 1 = a , x n + 1 = b x_1=a,x_{n+1}=b x1=a,xn+1=b,这样求定积分的问题就分解为下面的求和问题:

在这里插入图片描述

<2>数值积分的实现

  • 基于自适应辛普森方法:[l.n]=quad(filename, a, b, tol, trace)
  • 基于自适应Gauss-Lobatto方法:[l, n]=quadl(filename, a,b,tol, trace)

其中,filename是被积函数名;a和b分别是定积分的上限和下限,积分限[a, b]必须是有限的,不能是无穷大;tol用来控制积分精度,默认时取 t o l = 1 0 − 6 tol=10^{-6} tol=106;trace控制是否展现积分过程,若取非0则展现积分过程,取0则不展现,默认trace=0;返回参数I即定积分的值。n为被积函数的调用次数。

例子:

format long
f = @(x)4./(1+x.^2);
[I,n] = quad(f,0,1,1e-8)
[I,n] = quadl(f,0,1,1e-8)
(atan(1) - atan(0))*4

在这里插入图片描述

  • 基于全局自适应积分方法
    I=integral(filename,a,b)
    其中,I是计算得到的积分;filename是被积函数;a和b分别是定积分的下限和上限,积分限可以是无穷大
function f = f1(x)
f = 1./(x.*sqrt(1-log(x).^2));

>> I = integral(@f1,1,exp(1))

I =

   1.570796326795570

在这里插入图片描述

  • 基于自适应高斯-克朗罗德方法
    [I,err]=quadgk(filename, a, b)
    其中,err返回近似误差范围,其他参数的含义和用法与quad函数想相同。积分上下限可以是无穷大,也可以是复数。如果积分上下限为复数,则quadgk函数在复平面上求积分。

例子:

function f = f1(x)
f = sin(1./x)./x.^2;

I = quadgk(@f1,2/pi,+Inf)

I =

    1.0000

在这里插入图片描述

  • 基于梯形积分
    在这里插入图片描述

在这里插入图片描述

例子:

x = 1:6;
y = [6,8,11,7,5,2];
plot(x,y,'-ko');
grid on
axis([1,6,0,11]);
I1 = trapz(x,y)
%二者等价
I2 = sum(diff(x).*(y(1:end-1)+y(2:end))/2)

在这里插入图片描述
<3>多重定积分的数值求解
在这里插入图片描述

例子:

f2 = @(x,y)exp(-x.^2/2).*sin(x.^2+y);
I1 = quad2d(f2,-2,2,-1,1)
f3 = @(x,y,z)4*x.*z.*exp(-z.*z.*y-x.*x);
I2 = integral3(f3,0,pi,0,pi,0,1)

在这里插入图片描述

B 线性方程组求解

B.a 直接法

  • 高斯消去法
  • 列主元消去法
  • 矩阵的三角分解法

<1>利用左除运算符直接解法
Ax=b ====> x=A\b
注意:如果矩阵A是奇异的或者接近奇异的,则MATLAB会给出警告信息。

在这里插入图片描述

<2>利用矩阵分解求解线性方程组
1 LU分解
LU分解的基本方法:
在这里插入图片描述

MATLAB的LU分解函数:
[L,U]=lu(A):产生一个上三角阵U和一个变换形式的下三角阵L,使之满足A=LU。注意,这里的矩阵A必须是方阵。
[L,U,P]=lu(A):产生一个上三三角阵U和一个下三角阵L以及一个置换矩阵P,使之满足PA=LU。同样,矩阵A必须是方阵。

用LU分解求解线性方程组:
在这里插入图片描述

在这里插入图片描述

2 QR分解
3 Cholesky分解

B.b 迭代法

在这里插入图片描述

在这里插入图片描述

jacobi.m

%方程的解  迭代次数              迭代的初值 精度
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

在这里插入图片描述
gauseidel.m

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-L)\b;
 	y=B*x0+f;
 	n=1;
 	while norm(y-x0)>=ep
	  	x0=y;
 	 	y=B*x0+f;
  		n=n+1; 
 	end

例子:
在这里插入图片描述
有时候高斯-赛德尔迭代法求解线性方程组可能不收敛。
在这里插入图片描述
在这里插入图片描述

B.c 线性方程组应用举例

  • 平面桁架结构受力分析问题

  • 小行星运行轨道计算问题

    xi = [1.02,0.87,0.67,0.44,0.16];
    yi = [0.39,0.27,0.18,0.13,0.13];
    A = zeros(length(xi));
    for i = 1:length(xi)
        A(i,:) = [xi(i)*xi(i),2*xi(i)*yi(i),yi(i)*yi(i),2*xi(i),2*yi(i)];
    end
    b = -ones(length(xi),1);
    ai = A\b
    f = @(x,y)  2.4645*x.^2-0.8846*x.*y+6.4917*y.^2-1.3636*x-7.2016*y+1;
    h = ezplot(f,[-0.5,1.2,0,1.2]);
    

C 非线性方程求解与函数极值计算

C.a 非线性方程数值求解

<1>单变量非线性方程求解
函数的调用格式:x=fzero(filename, x0)
其中,filename是待求根方程左端的函数表达式,x0是初始值。

f = @(x) x-1./x+5;
x1 = fzero(f,-5)
x2 = fzero(f,1)
x3 = fzero(f,0.1)

在这里插入图片描述
<2>非线性方程组的求解
函数的调用格式为:
x=fsolve(filename, x0, option)
其中,x为返回的近似解,filename是待求根方程左端的函数表达式,x0是初值,option用于设置优化工具箱的优化参数,可以调用optimset函数来完成。 也可解例一

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)

在这里插入图片描述

C.b函数极值的计算

  • 极大值(最大值)
  • 极小值(最小值)

Matlab只考虑最小值问题的计算,如果要求 f ( x ) f(x) f(x)的最大值,可以求 − f ( x ) -f(x) f(x)的最小值
<1>无约束最优化问题
无约束最优化问题的一般描述为:
在这里插入图片描述
求最小值点和最小值的函数:
[xmin, fmin]=fminbnd(filename, x1, x2, option)
[xmin, fmin]=fminsearch(filename, x0, option)
[xmin, fmin]=fminunc(filename, x0, option)
fminbnd:一元函数
fminsearch:单纯形法。多元
fminunc:拟牛顿法。多元

其中,filename是定义的目标函数。第一个函数的输入变量x1、x2分别表示被研究区间的左右边界。后两个函数的输入变量x0是一个向量,表示极值点的初值。option为优化参数,可以通过optimset函数来设置。

f = @(x)x-1./x+5;
[xmin,fmin] = fminbnd(f,-10,-1)
[xmin,fmin] = fminbnd(f,1,10)

在这里插入图片描述

在这里插入图片描述

<2>有约束最优化问题
有约束最优化问题的一般描述为:
在这里插入图片描述
s . t . 是 s u b j e c t t o 的 缩 写 , 表 示 x 要 满 足 后 面 的 约 束 条 件 。 s.t.是subject to 的缩写,表示x要满足后面的约束条件。 s.t.subjecttox
约束条件可细化为:

  • 线性不等式约束
  • 线性等式约束
  • 非线性不等式约束
  • 非线性等式约束
  • x的上界和下界

求有约束条件下的最小值的函数为:
[xmin, fmin]=fmincon(filename, x0, A, b, Aeq, beq, Lbnd, Ubnd, Nonf, option)
其中,xmin, fmin,filename,x0和option的含义与求最小值函数相同。其余参数为约束条件,包括线性不等式约束、线性等式约束、x的下界和上界以及定义非线性约束的函数。如果某个约束不存在,则用空矩阵表示。

AX <= b(线性不等式约束,如果是约束是大于,则左右乘-1)
AeqX = beq(线性等式约束)
G(x) <= 0(非线性不等式约束)Nonf的两个输出值
Ceq(X) = 0(非线性等式约束)
Lbnd <= X <= Ubub(变量约束)

f = @(x)0.4*x(2)+x(1)^2+x(2)^2-x(1)*x(2)+1/30*x(1)^3;
x0 = [0.5;0.5];
A = [-1,-0.5;-0.5,-1];
b = [-0.4;-0.5];
lb = [0;0];
option = optimset('Display','off');
[xmin,fmin] = fmincon(f,x0,A,b,[],[],lb,[],[],option)

在这里插入图片描述

例题 建仓库送货问题

无约束最优化问题

a = [10,30,16.667,0.555,22.2221];
b = [10,50,29,29.888,49.988];
c = [10,18,20,14,25];
f = @(x)sum(c.*sqrt((x(1)-a).^2+(x(2)-b).^2));
[xmin,fmin] = fminsearch(f,[15,30])

有约束 y = x*x

%一般匿名函数返回值只有一个
function [c,h] = f1(x)
c= [];
h = x(2)-x(1)^2;

a = [10,30,16.667,0.555,22.2221];
b = [10,50,29,29.888,49.988];
c = [10,18,20,14,25];
f = @(x)sum(c.*sqrt((x(1)-a).^2+(x(2)-b).^2));
[xmin,fmin] = fmincon(f,[15,30],[],[],[],[],[],[],'f1')

D 常微分方程数值求解

D.a常微分方程数值求解的一般概念

在这里插入图片描述
在这里插入图片描述


D.b 常微分方程数值求解函数

在这里插入图片描述
在这里插入图片描述

MATLAB提供了多个求常微分方程初值问题数值解的函数,一般调用格式为:
[t,y]=solver(filename,tspan,y0,option)
其中,t和y分别给出时间向量和相应的数值解。
solver为求常微分方程数值解的函数。
filename是定义f(t,y)的函数名,该函数必须返回一个列向量。
tspan形式为[to,tf],表示求解区间。
y0是初始状态向量。
option是可选参数,用于设置求解属性,常用的属性包括相对误差值Re1To1(默认值是1 0 − 3 10^{-3}10
−3)和绝对误差值AbsTol(默认值是1 0 − 6 10^{-6}10 −6)。
常微分方程数值求解函数的统一命名格式:odennxx
ode是0rdinary Differential Equation的缩写,是常微分方程的意思。
nn是数字,代表所用方法的阶数。例如,ode23采用2阶龙格-库塔( R u n g e − K u t t a ) (Runge-Kutta)(Runge−Kutta)算法,用3阶公式做误差估计来调节步长,具有低等精度。
ode45采用4阶龙格-库塔算法,用5阶公式做误差估计来调节步长,具有中等精度。
xx是字母,用于标注方法的专门特征。
ode15s、ode23s中的“s”代表(Stiff),表示函数适用于刚性方程。
求常微分方程数值解的函数

例子:

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

在这里插入图片描述
例子:

f = @(t,x)[-2,0;0,1]*[x(2);x(1)];
[t,x] = ode45(f,[0,20],[1,0]);
subplot(2,2,1);
plot(t,x(:,2));
subplot(2,2,2);
plot(x(:,2),x(:,1));

在这里插入图片描述
在这里插入图片描述


D.c 刚性问题

在这里插入图片描述

在这里插入图片描述

lambda = 0.01;
f = @(t,y)y^2-y^3;
tic;[t,y] = ode45(f,[0,2/lambda],lambda);toc;
disp(['ode45计算的点数',num2str(length(t))]);

在这里插入图片描述

lambda = 1e-5;
f = @(t,y)y^2-y^3;
tic;[t,y] = ode15s(f,[0,2/lambda],lambda);toc;
disp(['ode45计算的点数',num2str(length(t))]);

在这里插入图片描述


F Lotka-Volterra模型

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-2vHBbzof-1625235977989)(C:\Users\myp\AppData\Roaming\Typora\typora-user-images\image-20210702153907321.png)]

第①问:兔子数量初始值x 1 ( 0 ) = 300,狐狸数量初始值x 2 ( 0 ) = 150

rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,30],[300,150])
subplot(1,2,1);
plot(t,x(:,1),'-',t,x(:,2),'-*');
legend('xl(t)','x2(t)');
xlabel('时间');
ylabel('物种数量');
grid on 
subplot(1,2,2);
plot(x(:,1),x(:,2))
grid on
1234567891011

6.61
我们可以看到,兔子和狐狸两者相互制约,当狐狸逐渐增多时,兔子逐渐减少,当狐狸增加到一定数量时,由于种群内部竞争激烈,又导致其数量减少,从而兔子天敌减少,进而导致兔子数量增加,当兔子增加到一定数量时,种间竞争加剧,从而又导致兔子数量减少,如此循环,相互制约发展。
**第②问:兔子初始值x 1 ( 0 ) = 15 ,狐狸初始值x 2 ( 0 ) = 22 **

rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,30],[15,22])
subplot(1,2,1);
plot(t,x(:,1),'-',t,x(:,2),'-*');
legend('xl(t)','x2(t)');
xlabel('时间');
ylabel('物种数量');
grid on 
subplot(1,2,2);
plot(x(:,1),x(:,2))
grid on

123456789101112

6.2

第③问:兔子初始值x 1 ( 0 ) = 102 ,狐狸初始值x 2 ( 0 ) = 198 。

rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,30],[102,198])
subplot(1,2,1);
plot(t,x(:,1),'-',t,x(:,2),'-*');
legend('xl(t)','x2(t)');
xlabel('时间');
ylabel('物种数量');
grid on 
subplot(1,2,2);
plot(x(:,1),x(:,2))
grid on
1234567891011

6.3

第④问:验证( 1 / λ , 2 / λ ) 是稳定平衡点。
取λ=0.01,所以稳定平衡点( 1 / λ , 2 / λ ) 即是( 100 , 200 ) ,以此点作为初值时,画出其图像。

rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,30],[100,200])
plot(t,x(:,1),'-o',t,x(:,2),'-*');
legend('x1(t)-兔子','x2(t)-狐狸');
xlabel('时间');
ylabel('物种数量');
123456

6.3
当将初始值变为( 98 , 195 ) (98,195)(98,195)时,即向下十分接近平衡点,画出其图像。

rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,30],[98,195])
plot(t,x(:,1),'-o',t,x(:,2),'-*');
legend('x1(t)-兔子','x2(t)-狐狸');
xlabel('时间');
ylabel('物种数量');
123456

6.4
两者图像没有交点,兔子和狐狸的数量分别以各自的平衡点进行周期性波动,但是两者的数量变化不是很大。
当将初始值变为( 70 , 150 ) (70,150)(70,150)时(向下偏离平衡点比较远时),画出其图像。
6.5

可以发现兔子和狐狸的数量变化比较剧烈,但是依然成周期性变化。
当将初始值变为( 900 , 1600 ) (900,1600)(900,1600)时(向上偏离平衡点十分远时),画出其图像。

rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,500],[900,1560])
plot(t,x(:,1),t,x(:,2));
legend('x1(t)-兔子','x2(t)-狐狸');
xlabel('时间');
ylabel('物种数量');
123456

6.6
可以发现,此时兔子和狐狸的数量变化相当的剧烈,但是依然存在周期。

F.a改进模型

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-997sU0h7-1625235977990)(C:\Users\myp\AppData\Roaming\Typora\typora-user-images\image-20210702153830684.png)]
第①问:在原模型下,狐狸数量和兔子数量的时间函数曲线

rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,50],[300,150])
plot(t,x(:,1),t,x(:,2));
legend('x1(t)-兔子','x2(t)-狐狸');
xlabel('时间');
ylabel('物种数量');
title('原模型下,狐狸和兔子数量的函数曲线');
1234567

6.7
第②问:在改进模型下,狐狸数量和兔子数量的时间函数曲线。

rabbitFox=@(t,x) [2*x(1)*(1-x(1)/400-0.005*x(2));...
    x(2)*(-1+0.01*x(1))];%...代表换行
[t,x]=ode45(rabbitFox,[0,50],[300,150]);
plot(t,x(:,1),t,x(:,2));
legend('x1(t)-兔子','x2(t)-狐狸');
xlabel('时间');
ylabel('物种数量');
title('改进模型下,狐狸和兔子数量的函数曲线');
12345678

在这里插入图片描述
第③问:在原模型下,狐狸数量相对兔子数量的关系曲线。

rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,30],[300,150])
plot(x(:,1),x(:,2));
xlabel('兔子数量');
ylabel('狐狸数量');
title('原模型下,狐狸数量相对于兔子数量的关系曲线');
123456

6.9
第④问:在改进模型下,狐狸数量相对于兔子数量的关系曲线

rabbitFox=@(t,x) [2*x(1)*(1-x(1)/400-0.005*x(2));...
    x(2)*(-1+0.01*x(1))];%...代表换行
[t,x]=ode45(rabbitFox,[0,50],[300,150]);
plot(x(:,1),x(:,2));
xlabel('时间');
ylabel('物种数量');
title('改进模型下,狐狸数量相对于兔子数量的关系曲线');
1234567

6.10

模型的进一步改进:
考虑狐狸的环境容纳量
考虑到兔子其他天敌和狐狸天敌的影响
考虑自然灾害、人为捕捉等因素

  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

追寻远方的人

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值