微分方程与欧拉法

微分方程概述

微分方程在各个领域应用颇多。

形如

的微分方程表示了系统的变化信息,

如果在加上初始条件(x0,y0),那么就可以求出系统整体随时间变化的信息。

可以说,正是微分方程将物理世界模型化

方向场与积分曲线

方向场(direction field)与积分曲线(integral curve)的关系,可以用下面的式子简要表示:

其中,当f(x,y),f′(x,y)在邻域内连续时,积分曲线不会相交也不会相切,解存在且唯一(exist and unique)。 下面,举函数y′=−x/y的方向场与积分曲线:

%下面的函数dirfield需要用到参考资料中的函数

%定义函数y'
f = @(x,y) -x/y

%方向场的一些简单可视化
ezplot(f,[-2,2,-2,2])
ezsurf(f,[-2,2,-2,2])
ezcontour(f,[-2,2,-2,2]); colorbar

%画出方向场与积分曲线
dirfield(f,-2:0.2:2,-2:0.2:2)
hold on 
for y0=-0.2:0.5:2 
  [ts,ys] = ode45(f,[-2,2],y0); plot(ts,ys) 
end 
title('dy/dx=-x/y的方向场与积分曲线')
hold off

微分方程的解析解法

微分方程的解析解法通常是将x,y分别移到等式的一边。 下面以

为例,移项后

所以有

进而有

最后解得:

其实,

就是根据微分方程y′=y在(0,1)(0,1)的初始条件下确定的。 使用matlab的解析解法为:

dsolve('Dy=2*y+1','x'
%输出为: (C2*exp(2*x))/2 - 1/2

%求解e^x
dsolve('Dy=y','y(0)=1','x')
%输出为: exp(x)

微分方程的数值解法

欧拉法

欧拉法的核心是,设定步长为h,然后已知y′和(x0,y0),根据下面方法迭代:

ODE数值解法的matlab程序为:

[xs,ys] = ode45(f,[-2,2],y0)

欧拉法的缺点

由上图可见,欧拉法存在一定的误差,并且误差会累计。当步长越小误差也就越小拟合效果越好。 这种情况下,误差和步长的关系是:

e∼c∗h

如果函数时而convex时而concave,这时候误差的变化便难以预测。

#---------------------------------------------------------
##凸函数
import numpy as np 
import scipy as sp
import matplotlib.pyplot as plt

#定义产生下一个点的函数
def nextPoint(x,y,h):
    xn = x + h
    slope = 2*x
    yn = y + h*slope
    return (xn,yn)
#定义产生点的生成器
def pointGenerator(x,y,h):
    while True:
        yield nextPoint(x,y,h)
        (x,y) = nextPoint(x,y,h)
#根据输入的起始终止点以及步长,输出可以用于画图的参数
def getXY(x,y,h):
    x1,y1=[],[]
    for i in pointGenerator(x,y,h):
        xi = i[0]
        yi = i[1]
        if xi>2.5:
            break
        else:
            x1.append(xi)
            y1.append(yi)
    x1.insert(0,-2)
    y1.insert(0,4)
    return (x1,y1)

#---------------------------------------------------------
#凹函数
#大部分和上面相同,只是将`nextPoint`函数重新定义
def nextPoint(x,y,h):
    xn = x + h
    slope = -2*x
    yn = y + h*slope
    return (xn,yn)

def pointGenerator(x,y,h):
    while True:
        yield nextPoint(x,y,h)
        (x,y) = nextPoint(x,y,h)

def getXY(x,y,h):
    x1,y1=[],[]
    for i in pointGenerator(x,y,h):
        xi = i[0]
        yi = i[1]
        if xi>2.5:
            break
        else:
            x1.append(xi)
            y1.append(yi)
    x1.insert(0,-2)
    y1.insert(0,-4)
    return (x1,y1)

x = np.arange(-2,2.1,0.1)
y = -x**2
x1,y1 = getXY(-2,-4,0.1)
x2,y2 = getXY(-2,-4,0.4)
x3,y3 = getXY(-2,-4,0.6)

plt.plot(x,y,'b--', linewidth=1,label='raw line')
plt.plot(x1,y1,'r',label='h=0.1')
plt.plot(x2,y2,'g',label='h=0.4')
plt.plot(x3,y3,'c',label='h=0.6')
plt.autoscale()
plt.xlim(-2.5,2.5)
plt.legend(loc='best')
plt.title('concave function with different h')
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

改进欧拉法之步长

步长的改进参考上文,步长越小误差越小。

改进欧拉法之斜率

核心是:计算斜率不只考虑当前的点,也考虑之后的点的斜率。 该方法一般被称作`runge-kutta`法,上文只用到一个斜率的被称为RK1,下面将要阐述的是RK2,同时在绝大多数数值计算工具中,`RK4`的使用最为广泛。

由上图可看,RK2的效果已经比RK1好太多的。

import numpy as np 
import scipy as sp
import matplotlib.pyplot as plt

#RK1
def nextPoint(x,y,h):
    xn = x + h
    slope = 2*x
    yn = y + h*slope
    return (xn,yn)


def pointGenerator(x,y,h):
    while True:
        yield nextPoint(x,y,h)
        (x,y) = nextPoint(x,y,h)


def getXY(x,y,h):
    x1,y1=[],[]
    for i in pointGenerator(x,y,h):
        xi = i[0]
        yi = i[1]
        if xi>2.5:
            break
        else:
            x1.append(xi)
            y1.append(yi)
    x1.insert(0,-2)
    y1.insert(0,4)
    return (x1,y1)


#RK2
def nextPoint2(x,y,h):
    xn = x + h
    slope = (2*x + 2*xn)/2
    yn = y + h*slope   
    return (xn,yn)


def pointGenerator2(x,y,h):
    while True:
        yield nextPoint2(x,y,h)
        (x,y) = nextPoint2(x,y,h)


def getXY2(x,y,h):
    x1,y1=[],[]
    for i in pointGenerator2(x,y,h):
        xi = i[0]
        yi = i[1]
        if xi>2.5:
            break
        else:
            x1.append(xi)
            y1.append(yi)
    x1.insert(0,-2)
    y1.insert(0,4)
    return (x1,y1)

#DATA
x = np.arange(-2,2.1,0.1)
y = x**2
x1,y1 = getXY(-2,4,1)
x2,y2 = getXY2(-2,4,1)
x3,y3 = getXY(-2,4,0.5)
x4,y4 = getXY2(-2,4,0.5)

#PLOT
plt.plot(x,y,'k', linewidth=1,label='raw line')
plt.plot(x4,y4,'r--',label='h=0.5 RK2')
plt.plot(x2,y2,'b--',label='h=1 RK2')
plt.plot(x3,y3,'c--',label='h=0.5 RK1')
plt.plot(x1,y1,'g--',label='h=1 RK1')

plt.autoscale()
plt.xlim(-2.5,2.5)
plt.legend(loc='best')
plt.title('convex function with different h and RKn')
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

参考资料

  1. MIT 18.03
  2. Using MATLAB for first order ODE

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏趣学算法

数据结构学习秘籍

网络上太多的同学吐槽被虐,如滔滔江水连绵不绝,数据结构太难了!真的很难吗?其实数据结构只是讲了三种:线性结构、树、图。到底难在哪里呢?通过调查了解大概有四个原因...

682
来自专栏C语言及其他语言

初学C语言的学习计划

背景:很多同学在学习C语言的过程中,常常会遇到这样的问题,即“教材看完了,知识点也懂,但写不出来程序”,这段时间,我们通过长期与有多年C语言研究经验的教授、教师...

3034
来自专栏一个会写诗的程序员的博客

计算机中的数学【阿贝尔-鲁菲尼定理】五次方程的根

这个定理以保罗·鲁菲尼和尼尔斯·阿贝尔命名。前者在1799年给出了一个不完整的证明,后者则在1824年给出了完整的证明。埃瓦里斯特·伽罗瓦创造了群论,独立地给出...

732
来自专栏鹅厂优文

游戏人工智能 读书笔记 (四) AI算法简介——Ad-Hoc 行为编程

本书英文版: Artificial Intelligence and Games - A Springer Textbook

15610
来自专栏tkokof 的技术,小趣及杂念

对于"矩阵连乘问题"的一点想法

在算法设计的学习中,每到“动态规划”一节,一般都会涉及到“矩阵连乘”问题(例如《Algorithms》,中文译名《算法概论》),可想而知该题的经典程度 :)

613
来自专栏数据结构与算法

19:救援

19:救援 总时间限制: 1000ms 内存限制: 65536kB描述 救生船从大本营出发,营救若干屋顶上的人回到大本营,屋顶数目以及每个屋顶的坐标  和人数...

3399
来自专栏专知

【LeetCode 204】关关的刷题日记40 Number of Boomerangs

关关的刷题日记40 – Leetcode 447. Number of Boomerangs 题目 Given n points in the plane th...

3244
来自专栏IT派

Python 实现简单的导弹自动追踪

自动追踪算法,在我们设计2D射击类游戏时经常会用到,这个听起来很高大上的东西,其实也并不是军事学的专利,在数学上解决的话需要去解微分方程,

843
来自专栏机器学习原理

我的机器学习微积分篇观点函数从极限到导数导数的应用偏导数从方向导数到梯度

前言: 没想到还能在此生再次用到大学习学习的高数,线性代数和概率论,如果上天给我再来一次的机会,我一定往死了学习这三门课。 观点 与机器学习相关的微积分...

3765
来自专栏懒人开发

(2.1)James Stewart Calculus 5th Edition:The Tangent and Velocity Problems

这里提到 Tangent 起源于 拉丁文, 意思是 touching 也就是曲线对应点位置当前的方向

842

扫码关注云+社区