一、欧拉法

欧拉法：是一种微分方程的数值计算方法，其基本思想是迭代。所谓迭代，就是逐次替代，最后求出所要求的解，并达到一定的精度。误差可以很容易地计算出来。下面分别介绍欧拉法求解一阶微分方程，二阶微分方程的matlab程序实现。

1.一阶微分方程

{dy(t)dt=f(y(t),t)y(0)=y0\begin{cases} {dy(t) \over dt} =f(y(t),t) \\ y(0)=y_0 \end{cases} {dtdy(t)=f(y(t),t)y(0)=y0
用欧拉法实现的步骤如下：

第一步：在给定的区间[0，T]进行分割，将区间T进行N等分（Δt），分别记为t1,t2,t3,…tn+1；

第二步：根据给定条件，给定一阶导数(导数代表自变量的变化率)和初始值y0，利用talyor展开可以把第n个点t的值y(tn)用f’(tn-1)和tn-1表示出来；
y(tn+1)=y(tn)+f(y(t),t)Δt+ο(Δt2)y(t_n+_1) =y(t_n)+f(y(t),t)\Delta t+\omicron(\Delta t^2) y(tn+1)=y(tn)+f(y(t),t)Δt+ο(Δt2)

第三步：近似，忽略掉高阶无穷小（拉格朗日余项）进行近似，根据一致的y0进行不断迭代；
y(tn+1)=y(tn)+f(y(t),t)Δty(t_n+_1) =y(t_n)+f(y(t),t)\Delta t y(tn+1)=y(tn)+f(y(t),t)Δt
Matlab实现如下：
以{dy(t)dt=ety(0)=1\begin{cases} {dy(t) \over dt} =e^t \\ y(0)=1 \end{cases}{dtdy(t)=ety(0)=1为例

clear
y0=1;
T=2;
dt=T/N;%分割为1000份（或更多），用于控制精度
t=0:dt:2;
y1=zeros(1,length(t)-1);
y1(1)=y0;
for n=1:length(t)-1
y1(n+1)=y1(n)+exp(t(n))*dt;
end
plot(t,y1,t,exp(t))
可以看出和真实曲线基本完全重合，读者可以利用其他已知函数来验证。

2.二阶微分方程

{d2y(t)dt2=f(y(t),t,y′(t))y(0)=y0dy(0)dt=v0\begin{cases} {d^2y(t) \over dt^2} =f(y(t),t,y'(t)) \\ y(0)=y_0 \\ {dy(0) \over dt} =v_0 \end{cases} ⎩⎪⎨⎪⎧dt2d2y(t)=f(y(t),t,y′(t))y(0)=y0dtdy(0)=v0
具体的原理可以参考‘一阶微分方程’的做法，本系列文章主要用于物理问题的模拟，所以接下来用简谐运动为例来说明，
{d2y(t)dt2=−kmy(t)y(0)=y0dy(0)dt=v0\begin{cases} {d^2y(t) \over dt^2} =-\text{\(\frac k m\)} y(t)\\ y(0)=y_0 \\ {dy(0) \over dt} =v_0 \end{cases} ⎩⎪⎨⎪⎧dt2d2y(t)=−mky(t)y(0)=y0dtdy(0)=v0
这个微分方程就有了物理意义，即给定初始位置和速度条件下的简谐运动，可以根据提供的参量通过欧拉法把整个运动定下来，得到任意时刻的速度与位移！

clear
k=1;
m=1;
y0=0;
v0=1;
T=20;
N=100000;%用于控制精度,N越大越接近真实结果
dt=T/N;
t=0:dt:T;
y1=zeros(1,length(t));%产生位移向量
y2=zeros(1,length(t));%产生速度向量
y1(1)=y0;
y2(1)=v0;
for n=1:length(t)-1
y1(n+1)=y1(n)+y2(n)*dt;
y2(n+1)=y2(n)+(-k/m)*y1(n)*dt;
end
plot(t,y1,t,y2)
蓝色线代表y1-t图（位移时间），橙色代表y2-t图（速度时间），读者可以调试k，m，y0，v0得到不同的简谐运动。

相空间

尚若在上边程序中最终执行：

plot(y1,y2)
axis equal%用于使两个坐标轴等刻度
%获得以下图像
通常我们是以x，y为横纵坐标，但也可以用x（位移）和p（动量）分别作为横纵坐标，得到以上图像，即为在相空间中的图像，但是为什么是个正圆呢？
这是因为能量守恒以及k，m的取值决定：
首先在这个例子中没有阻尼，所以能量守恒，即振子的动能和势能合为一个常数，1/2mv^2 + 1/2kx^2=C(常数)；
又因为k=m=1，所以上边方程可以化成圆的方程的一般形式所以最终为正圆；
若改变k，m的取值，可想而之最终的结果是个椭圆，读者可以自行尝试。

总结：未完待续！

数值方法求解微分方程相关推荐

栅格单元值的选取方法_计算机求解微分方程的六大数值计算方法
今天介绍关于计算机求解微分方程的六大数值计算的方法. 1.有限元法有限元方法的基础是变分原理和加权余量法,其基本求解思想是把计算域划分为有限个互不重叠的单元,在每个单元内,选择一些合适的节点作为求解 ...
C#，数值计算，求解微分方程的吉尔（Gear）四阶方法与源代码
微分方程微分方程,是指含有未知函数及其导数的关系式.解微分方程就是找出未知函数. 微分方程是伴随着微积分学一起发展起来的.微积分学的奠基人Newton和Leibniz的著作中都处理过与微分方程有关的 ...
matlab:使用欧拉方法求解微分方程
书目:<数值方法(matlab版)>,作者周璐等译 %欧拉方法求解微分方程function [t,y] = my_euler(f, t0, tf , y0, h)%f-函数: t0,tf: ...
AI攻破高数核心，1秒内精确求解微分方程、不定积分，性能远超Matlab
栗子鱼羊发自海边边量子位报道 | 公众号 QbitAI 大家都知道,AI (神经网络) 连加减法这样的简单算术都做不好: 可现在,AI已经懂得微积分,把魔爪伸向你最爱的高数了. 它不光会求 ...
MATLAB从入门到精通-欧拉法与梯形法求解微分方程（含MATLAB源码）
前言以下是我为大家准备的几个精品专栏,喜欢的小伙伴可自行订阅,你的支持就是我不断更新的动力哟! MATLAB-30天带你从入门到精通 MATLAB深入理解高级教程(附源码) tableau可视化数据 ...
matlab求微分方程同届,Matlab学习——求解微分方程(组)
介绍: 1.在 Matlab 中,用大写字母 D 表示导数,Dy 表示 y 关于自变量的一阶导数,D2y 表示 y 关于自变量的二阶导数,依此类推.函数 dsolve 用来解决常微分方程(组)的求解问 ...
AI攻破高数核心，1秒内精确求解微分方程、不定积分
大家都知道,AI (神经网络) 连加减法这样的简单算术都做不好: 可现在,AI已经懂得微积分,把魔爪伸向你最爱的高数了. 它不光会求不定积分: 还能解常微分方程: 一阶二阶都可以. 这是Faceboo ...
python求解微分方程_python能解微分方程吗
Sympy Sympy是一个数学符号库(sym代表了symbol,符号),包括了积分,微分方程等各种数学运算方法,为python提供了强大的数学运算支持.对于图像来说,虽然都是做离散的计算,操作最多的 ...
龙格库塔法和欧拉法求解微分方程的比较
文章目录计算机如何理解连续系统的动态特性? 欧拉法求解微分方程龙格库塔法求解微分方程 MATLAB代码编写和仿真效果计算机如何理解连续系统的动态特性? 一般连续系统的动态特性可以由一个微分方程, ...
matlab求解微分方程6,牛津大学出版社数学应用软件作业6 用Matlab求解微分方程(组)的解析解和数值解:经典爱情语录大全...
漳州理工职业学院-酒会礼仪注:上机作业文件夹以自己的班级姓名学号命名,文件夹包括如下上机报告和Matlab程序. 上机报告模板如下: 佛山科学技术学院上机报告课程名称数学应用软件上机项 ...

数值方法求解微分方程

文章目录

一、欧拉法

1.一阶微分方程

2.二阶微分方程

相空间

总结：未完待续！

数值方法求解微分方程相关推荐

最新文章

热门文章