您好,欢迎访问代理记账网站
移动应用 微信公众号 联系我们

咨询热线 -

电话 15988168888

联系客服
  • 价格透明
  • 信息保密
  • 进度掌控
  • 售后无忧

matlab求解微分方程的数值解

简 介:前面介绍了微分方程的解析解方法,同时也指出很多非线性微分方程是不存在解析解法的,需要使用数值解法对之进行研究。本文着重讨论基于 MATLAB/Simulink语言的各类微分方程的数值解方法。

关键词 微分方程、数值解、MATLAB

§01


一般微分方程的数值解法很大一类是关于微分方程初值问题的数值解法,这类问题需要用一阶显式的微分方程组来描述为

x ˙ ( t ) = f ( t , x ( t ) ) \dot{\boldsymbol{x}}(t)=\boldsymbol{f}(t, \boldsymbol{x}(t)) x˙(t)=f(t,x(t))

其中, x T ( t ) = [ x 1 ( t ) , x 2 ( t ) , ⋯   , x n ( t ) ] \boldsymbol{x}^{\mathrm{T}}(t)=\left[x_{1}(t), x_{2}(t), \cdots, x_{n}(t)\right] xT(t)=[x1(t),x2(t),,xn(t)]称为状态向量, f T ( ⋅ ) = [ f 1 ( ⋅ ) , f 2 ( ⋅ ) , ⋯   , f n ( ⋅ ) ] \boldsymbol{f}^{\mathrm{T}}(\cdot)=\left[f_{1}(\cdot), f_{2}(\cdot), \cdots, f_{n}(\cdot)\right] fT()=[f1(),f2(),,fn()]可以是任意非线性函数。所谓初值问题是指,若已知初始状态 x 0 = [ x 1 ( 0 ) , ⋯   , x n ( 0 ) ] T \boldsymbol{x}_{0}=\left[x_{1}(0), \cdots, x_{n}(0)\right]^{\mathrm{T}} x0=[x1(0),,xn(0)]T,用数值求解方法求出在某个时间区间 t ∈ [ 0 , t f ] t \in\left[0, t_{\mathrm{f}}\right] t[0,tf]内各个时刻状态变量 x ( t ) \boldsymbol{x}(t) x(t) 的数值,这里 t f t_{\mathrm{f}} tf 又称为终止时间。

其他类型微分方程求解必须先转换:

1. 单个高阶常微分方程处理方法

一个高阶常微分方程的一般形式为:

y ( n ) = f ( t , y , y ′ , ⋯   , y ( n − 1 ) ) y^{(n)}=f\left(t, y, y^{\prime}, \cdots, y^{(n-1)}\right) y(n)=f(t,y,y,,y(n1))

输出变量的各阶导数初始值为:

y ( 0 ) ,    y ′ ( 0 ) ,    ⋯   ,    y ( n − 1 ) ( 0 ) y(0), ~~y^{\prime}(0), ~~\cdots,~~ y^{(n-1)}(0) y(0),  y(0),  ,  y(n1)(0)

选择一组状态变量:

x 1 = y ,    x 2 = y ′ ,    ⋯   ,    x n = y ( n − 1 ) x_{1}=y, ~~x_{2}=y^{\prime}, ~~\cdots,~~ x_{n}=y^{(n-1)} x1=y,  x2=y,  ,  xn=y(n1)

原高阶常微分方程模型变换为一阶微分方程组:

{ x 1 ′ = x 2 x 2 ′ = x 3 ⋮ x n ′ = f ( t , x 1 , x 2 , ⋯   , x n ) \left\{\begin{aligned} x_{1}^{\prime} &=x_{2} \\ x_{2}^{\prime} &=x_{3} \\ & \vdots \\ x_{n}^{\prime} &=f\left(t, x_{1}, x_{2}, \cdots, x_{n}\right) \end{aligned}\right. x1x2xn=x2=x3=f(t,x1,x2,,xn)

初值为:

x 1 ( 0 ) = y ( 0 ) ,    x 2 ( 0 ) = y ′ ( 0 ) ,    ⋯   ,    x n ( 0 ) = y ( n − 1 ) ( 0 ) x_{1}(0)=y(0), ~~x_{2}(0)=y^{\prime}(0), ~~\cdots, ~~x_{n}(0)=y^{(n-1)}(0) x1(0)=y(0),  x2(0)=y(0),  ,  xn(0)=y(n1)(0)

2. 高阶常微分方程组的变换方法

多元高阶常微分方程组的处理

{ x ( m ) = f ( t , x , x ′ , ⋯   , x ( m − 1 ) , y , ⋯   , y ( n − 1 ) ) y ( n ) = g ( t , x , x ′ , ⋯   , x ( m − 1 ) , y , ⋯   , y ( n − 1 ) ) \left\{\begin{array}{l}x^{(m)}=f\left(t, x, x^{\prime}, \cdots, x^{(m-1)}, y, \cdots, y^{(n-1)}\right) \\ \\ y^{(n)}=g\left(t, x, x^{\prime}, \cdots, x^{(m-1)}, y, \cdots, y^{(n-1)}\right)\end{array}\right. x(m)=f(t,x,x,,x(m1),y,,y(n1))y(n)=g(t,x,x,,x(m1),y,,y(n1))

状态变量的选择不唯一,建议:选择如下状态变量

x 1 = x ,   x 2 = x ′ ,   ⋯   ,   x m = x ( m − 1 ) ,   x m + 1 = y ,   x m + 2 = y ′ ,   ⋯   ,   x m + n = y ( n − 1 ) x_{1}=x, ~x_{2}=x^{\prime}, ~\cdots, ~x_{m}=x^{(m-1)}, \\ \ \\ x_{m+1}=y, ~x_{m+2}=y^{\prime},~ \cdots, ~x_{m+n}=y^{(n-1)} x1=x, x2=x, , xm=x(m1), xm+1=y, xm+2=y, , xm+n=y(n1)

新的状态方程

{ x 1 ′ = x 2 ⋮ x m ′ = f ( t , x 1 , x 2 , ⋯   , x m + n ) x m + 1 ′ = x m + 2 ⋮ x m + n ′ = g ( t , x 1 , x 2 , ⋯   , x m + n ) \left\{\begin{aligned} x_{1}^{\prime}=& x_{2} \\ & \vdots \\ x_{m}^{\prime}=& f\left(t, x_{1}, x_{2}, \cdots, x_{m+n}\right) \\ x_{m+1}^{\prime} &=x_{m+2} \\ & \vdots \\ x_{m+n}^{\prime} &=g\left(t, x_{1}, x_{2}, \cdots, x_{m+n}\right) \end{aligned}\right. x1=xm=xm+1xm+nx2f(t,x1,x2,,xm+n)=xm+2=g(t,x1,x2,,xm+n)

可以描述该方程,然后用 ode45 等求解。

§02 分方程求解的误差与步长问题


1. 不能无限制地减小步长 h h h的值的两条原因

  • 减慢计算速度
  • 增加累积误差

2. 在对微分方程求解过程中应采取的三个措施

  • 选择适当的步长
  • 改进近似算法精度
  • 采用变步长方法

§03 数调用格式(ode45)


[t, x] = ode45(fun,[t0, tf], x0)  % 直接求解
[t, x] = ode45(fun,[t0, tf], x0, options)  % 带有控制选项
[t, x] = ode45(fun,[t0, tf], x0, options, p1, p2, ...)  % 带有附加参数

当然,还存在其他求解函数,如ode15sode23等,不同的算法(函数)适合于不同类型的微分方程。

§04 分方程的MATLAB描述


描述需要求解的微分方程组

   f u n c t i o n     x d = f u n ( t , x ) \rm{function}~~~ \boldsymbol{x}_{d}= fun (t, \boldsymbol{x}) function   xd=fun(t,x)

   f u n c t i o n     x d = f u n ( t , x , p 1 , p 2 , ⋯   ) \rm{function} ~~~ \boldsymbol{x}_{d}= fun \left(t, \mathrm{x}, p_{1}, p_{2}, \cdots\right) function   xd=fun(t,x,p1,p2,)

修改控制变量

options = odeset('RelTol', 1e-7) ;
options = odeset;  options.RelTol = 1e-7;

§05 用举例:Lorenz方程


例1:无参数

题目: 求解下列Lorenz模型

{ x ˙ 1 ( t ) = − β x 1 ( t ) + x 2 ( t ) x 3 ( t ) x ˙ 2 ( t ) = − ρ x 2 ( t ) + ρ x 3 ( t ) x ˙ 3 ( t ) = − x 1 ( t ) x 2 ( t ) + σ x 2 ( t ) − x 3 ( t ) \left\{\begin{array}{l}\dot{x}_{1}(t)=-\beta x_{1}(t)+x_{2}(t) x_{3}(t) \\ \\ \dot{x}_{2}(t)=-\rho x_{2}(t)+\rho x_{3}(t) \\ \\ \dot{x}_{3}(t)=-x_{1}(t) x_{2}(t)+\sigma x_{2}(t)-x_{3}(t)\end{array}\right. x˙1(t)=βx1(t)+x2(t)x3(t)x˙2(t)=ρx2(t)+ρx3(t)x˙3(t)=x1(t)x2(t)+σx2(t)x3(t)

式中参数为 β = 8 / 3 ,   ρ = 10 ,   σ = 28 \beta=8 / 3, ~\rho=10, ~\sigma=28 β=8/3, ρ=10, σ=28

初始条件为 x 1 ( 0 ) = x 2 ( 0 ) = 0 ,    x 3 ( 0 ) = ϵ ,    i . e . ,   ϵ = 1 0 − 10 x_{1}(0)=x_{2}(0)=0, ~~x_{3}(0)=\epsilon,~~ \rm{i.e.}, ~\epsilon=10^{-10} x1(0)=x2(0)=0,  x3(0)=ϵ,  i.e., ϵ=1010

求解:

1. 写出标准型

x ′ ( t ) = [ − β x 1 ( t ) + x 2 ( t ) x 3 ( t ) − ρ x 2 ( t ) + ρ x 3 ( t ) − x 1 ( t ) x 2 ( t ) + σ x 2 ( t ) − x 3 ( t ) ] , x ( 0 ) = [ 0 0 ϵ ] \boldsymbol{x}^{\prime}(t)=\left[\begin{array}{c}-\beta x_{1}(t)+x_{2}(t) x_{3}(t) \\ \\ -\rho x_{2}(t)+\rho x_{3}(t) \\ \\ -x_{1}(t) x_{2}(t)+\sigma x_{2}(t)-x_{3}(t)\end{array}\right], \quad \boldsymbol{x}(0)=\left[\begin{array}{l}0 \\ \\ 0 \\ \\ \epsilon\end{array}\right] x(t)=βx1(t)+x2(t)x3(t)ρx2(t)+ρx3(t)x1(t)x2(t)+σx2(t)x3(t),x(0)=00ϵ

2. 微分方程的MATLAB描述

  • M-文件描述
function y = lorenzeq(t, x)
	y = [-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3); -x(1)*x(2)+28*x(2)-x(3)];
  • 匿名函数
f = @(t,x)[-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3); -x(1)*x(2)+28*x(2)-x(3)];

3. 微分方程求解

  • 匿名函数
t_final = 100; 
x0 = [0; 0; 1e-10]; 
f = @(t,x)[-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3); -x(1)*x(2)+28*x(2)-x(3)];
[t,x] = ode45(f, [0,t_final], x0); 
plot(t,x), figure; plot3(x(:,1),x(:,2),x(:,3));
  • M-文件描述
t_final = 100; 
x0 = [0; 0; 1e-10]; 
f = @(t,x)[-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3); -x(1)*x(2)+28*x(2)-x(3)];
[t,x] = ode45(@lorenzeq, [0,t_final], x0); 
plot(t,x), figure; plot3(x(:,1),x(:,2),x(:,3));

function y = lorenzeq(t, x)
	y = [-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3); -x(1)*x(2)+28*x(2)-x(3)];
end
▲ 图 1

▲ 图 2

例2:带有附加参数

引入附加参数的目的: 如果参数 β \beta β, ρ \rho ρ, σ \sigma σ 改变,不需要修改原函数。

重新求解Lorenz方程

方式一

f = @(t,x,beta,rho,sigma)[-beta*x(1) + x(2)*x(3); -rho*x(2) + rho*x(3); -x(1)*x(2) + sigma*x(2) - x(3)];
t_final=100; 
x0 = [0; 0; 1e-10]; 
b1 = 8/3; r1 = 10; s1 = 28; 
[t,x]=ode45(f, [0,t_final], x0, [], b1, r1, s1);

方式二

beta = 2; rho = 5; sigma = 20; 
f = @(t,x)[-beta*x(1) + x(2)*x(3); -rho*x(2) + rho*x(3); -x(1)*x(2) + sigma*x(2) - x(3)];
t_final = 100;
x0 = [0; 0; 1e-10]; 
[t2,x2] = ode45(f, [0,t_final], x0);

分享:

低价透明

统一报价,无隐形消费

金牌服务

一对一专属顾问7*24小时金牌服务

信息保密

个人信息安全有保障

售后无忧

服务出问题客服经理全程跟进