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))
其中, xT(t)=[x1(t),x2(t),⋯,xn(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)]称为状态向量,fT(⋅)=[f1(⋅),f2(⋅),⋯,fn(⋅)]\boldsymbol{f}^{\mathrm{T}}(\cdot)=\left[f_{1}(\cdot), f_{2}(\cdot), \cdots, f_{n}(\cdot)\right]fT(⋅)=[f1(⋅),f2(⋅),⋯,fn(⋅)]可以是任意非线性函数。所谓初值问题是指,若已知初始状态x0=[x1(0),⋯,xn(0)]T\boldsymbol{x}_{0}=\left[x_{1}(0), \cdots, x_{n}(0)\right]^{\mathrm{T}}x0=[x1(0),⋯,xn(0)]T,用数值求解方法求出在某个时间区间t∈[0,tf]t \in\left[0, t_{\mathrm{f}}\right]t∈[0,tf]内各个时刻状态变量 x(t)\boldsymbol{x}(t)x(t) 的数值,这里 tft_{\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(n−1))
输出变量的各阶导数初始值为:
y(0),y′(0),⋯,y(n−1)(0)y(0), ~~y^{\prime}(0), ~~\cdots,~~ y^{(n-1)}(0)y(0), y′(0), ⋯, y(n−1)(0)
选择一组状态变量:
x1=y,x2=y′,⋯,xn=y(n−1)x_{1}=y, ~~x_{2}=y^{\prime}, ~~\cdots,~~ x_{n}=y^{(n-1)}x1=y, x2=y′, ⋯, xn=y(n−1)
原高阶常微分方程模型变换为一阶微分方程组:
{x1′=x2x2′=x3⋮xn′=f(t,x1,x2,⋯,xn)\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.⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧x1′x2′xn′=x2=x3⋮=f(t,x1,x2,⋯,xn)
初值为:
x1(0)=y(0),x2(0)=y′(0),⋯,xn(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(n−1)(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(m−1),y,⋯,y(n−1))y(n)=g(t,x,x′,⋯,x(m−1),y,⋯,y(n−1))
状态变量的选择不唯一,建议:选择如下状态变量
x1=x,x2=x′,⋯,xm=x(m−1),xm+1=y,xm+2=y′,⋯,xm+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(m−1), xm+1=y, xm+2=y′, ⋯, xm+n=y(n−1)
新的状态方程
{x1′=x2⋮xm′=f(t,x1,x2,⋯,xm+n)xm+1′=xm+2⋮xm+n′=g(t,x1,x2,⋯,xm+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+1′xm+n′x2⋮f(t,x1,x2,⋯,xm+n)=xm+2⋮=g(t,x1,x2,⋯,xm+n)
可以描述该方程,然后用 ode45
等求解。
§02 微分方程求解的误差与步长问题
1. 不能无限制地减小步长hhh的值的两条原因
- 减慢计算速度
- 增加累积误差
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, ...) % 带有附加参数
当然,还存在其他求解函数,如ode15s
、ode23
等,不同的算法(函数)适合于不同类型的微分方程。
§04 微分方程的MATLAB描述
描述需要求解的微分方程组
functionxd=fun(t,x)\rm{function}~~~ \boldsymbol{x}_{d}= fun (t, \boldsymbol{x})function xd=fun(t,x)
functionxd=fun(t,x,p1,p2,⋯)\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)=−βx1(t)+x2(t)x3(t)x˙2(t)=−ρx2(t)+ρx3(t)x˙3(t)=−x1(t)x2(t)+σx2(t)−x3(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
初始条件为 x1(0)=x2(0)=0,x3(0)=ϵ,i.e.,ϵ=10−10x_{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., ϵ=10−10
求解:
1. 写出标准型
x′(t)=[−βx1(t)+x2(t)x3(t)−ρx2(t)+ρx3(t)−x1(t)x2(t)+σx2(t)−x3(t)],x(0)=[00ϵ]\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);
matlab求解微分方程的数值解相关推荐
- matlab求解微分方程6,牛津大学出版社数学应用软件作业6 用Matlab求解微分方程(组)的解析解和数值解:经典爱情语录大全...
漳州理工职业学院-酒会礼仪 注:上机作业文件夹以自己的班级姓名学号命名,文件夹包括如下上机报告和Matlab程序. 上机报告模板如下: 佛山科学技术学院 上 机 报 告 课程名称 数学应用软件 上机项 ...
- matlab解无解析解微分方程组,数学应用软件作业6 用Matlab求解微分方程(组)的解析解和数值解...
<数学应用软件作业6 用Matlab求解微分方程(组)的解析解和数值解>由会员分享,可在线阅读,更多相关<数学应用软件作业6 用Matlab求解微分方程(组)的解析解和数值解(12页 ...
- matlab求解全微分函数,利用MATLAB求解微分方程的方法探索
引言 科学问题和工程问题经常需要求取微分方程的解,MATLAB 的强大数值运算和符号运算能力,能够方便地进行各种解析运算,是方便实用.功能强大的数学软件之一. 1线性微分方程求解 1.1线性常微分方程 ...
- Matlab求解微分方程数值解
有三种方法求解微分方程数值解: 欧拉法 改进欧拉法 龙格库塔法 接下来用一个练习来对比这三种求解方法. 问题描述:用改进的Euler方法.MATLAB的ode45命令分别求下列初值问题的数值解,并画图 ...
- Matlab求微分方程的数值解
注:首先计算微分方程的解析解,如果发现没有解析解,再用数值解 一.Matlab中求微分方程的数值解函数 [x,y]=solver('f',ts,x0,options) 1)x代表自变量 2)y代表函数 ...
- matlab求解微分方程解析解
解析解 1求解微分方程方程组用函数"dsolve" dsolve('方程1','方程2'-'方程n','初始条件','自变量') clc,clear %% %求解du/dt=i+u ...
- 数学建模学习(28):又一夜没睡,爆肝整理所有类型matlab求解微分方程+案例实战,学不会来砍我
文章目录 前言 dsolve语法 参数详解 求解一阶微分方程 求解二阶微分方程 求解有条件的微分方程 将输出分配给函数或变量 求微分方程的显式和隐式解 当未找到显式解决方案时查找隐式解决方案 求微分方 ...
- Matlab答疑五:使用微分定义求解微分方程的数值解
1.题目 解微分方程 dydt=sin(y)+t,其中t=0时y=0,并绘图. 说明,一般对dydt的求解方法为:y(t+dt)=y(t)+dydt(t)*dt 2.方法 除了题目给出方法:使用定义求 ...
- matlab偏导数方程,[转载]Matlab求解微分方程(2)——偏微分方程的求解
从写完上一篇常微分方程的求解到现在已经很长时间了,这周也一直忙于报到的各种事宜,无暇坐下来写些东西,趁着这个周末,终于完成了这个姊妹篇. 对于偏微分方程的求解,Matlab提供了两种工具.第一种是pd ...
最新文章
- Linux目录是否是否为空,在Linux上使用C检查目录是否为空
- SPOJ 375. Query on a tree (树链剖分)
- [转发]项目修复-把有麻烦的项目带向成功
- 【进程通信】Socket
- flutter怎么手动刷新_flutter局部刷新的实现示例
- Android 创建,删除,检测桌面快捷方式
- kafka吞吐量高的原因
- mysql 逗号金额比较,如何使用MySQL比较两个逗号分隔的字符串列表
- dataframe 排序_疯狂Spark之DataFrame创建方式详解一(九)
- 微服务_SpringCloud微服务架构实战:高并发微服务架构设计
- WebService开发常用功能详解
- win7 IIS服务启动和停止
- Atitit 架构师之道 attilax著 1.1. 认和评估系统需求,	2 1.2. 给出开发规范	2 1.3. ,搭建系统实现的核心构架,	2 1.4. 扫清主要难点的技术人员	2 1.5. 核
- 宝软网java软件下载_手机游戏怎么下载
- 怎么用dos系统进入服务器,怎么用DOS命令方式启动系统服务
- c语言中begin用法,C++ set cbegin() 使用方法及示例
- 斑马打印机打印中文乱码的问题
- 石墨笔记,熊掌记和 Effie 哪个更适合采编?
- 换IP软件能否实现定时切换IP?一起来验证
- matlab实现图像DCT变换