# 方向场与积分曲线

```%下面的函数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```

# 微分方程的解析解法

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

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

# 微分方程的数值解法

## 欧拉法

ODE数值解法的matlab程序为：

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

## 欧拉法的缺点

e∼c∗h

```#---------------------------------------------------------
##凸函数
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()```

## 改进欧拉法之斜率

```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()```

# 参考资料

99 篇文章32 人订阅

0 条评论

## 相关文章

8700

31610

1392

### 2015美团校招部分笔试题

http://blog.csdn.net/lanxuezaipiao/article/details/41774539

1042

40111

701

1214

28510

### tf28: 手写汉字识别

MNIST手写数字数据集通常做为深度学习的练习数据集，这个数据集恐怕早已经被大家玩坏了。本帖就介绍一个和MNIST类似，同时又适合国人练习的数据集-手写汉字数...

4419

### scikit-learn 梯度提升树(GBDT)调参小结

在梯度提升树(GBDT)原理小结中，我们对GBDT的原理做了总结，本文我们就从scikit-learn里GBDT的类库使用方法作一个总结，主要会关注调参...

1243