2021-11-12

解常微分方程

1. 利用sympy的odient解ode
   以二阶为例
   $$f(x)=\{\begin{array}{rl}{y}^{\prime \prime }-y+x& =0\\ y(0)=y(1)& =0\end{array}$$
   已知答案：
   $$y^{\prime }=x-\frac{sinh(x)}{sinh(1)}$$

from scipy.integrate import odeint
from sympy.abc import t
import numpy as np
import matplotlib.pyplot as plt
import matplotlib



#定义一个方程组（微分方程组）
def pfun(y,x):
	y1,y2=y; #让'y'成为一个[y1',y2']的向量 所以将等式左边都化为1阶微分是很重要的
	return np.array([y2,y1-x]) #返回的是等式右边的值
x=np.arange(0,1,0.01) #创建自变量序列
soli=odeint(pfun,[0.0,0.149],x) #求数值解
plt.rc('font',size=16); plt.rc('font',family='SimHei')
plt.plot(x,soli[:,0],'r*',label="数值解")
plt.rcParams['axes.unicode_minus']=False
plt.plot(x,x-np.sinh(x)/np.sinh(1),'g',label="符号解曲线")
plt.legend()
plt.show()

针对于odeint必须有初始值，方程组未给一阶导初始值，这里利用答案反推`

import numpy as np
import sympy as s


def f(x):
    return x-np.sinh(x)/np.sinh(1)#x-np.sinh(x)/np.sinh(1)
x = s.symbols("x")
print(s.diff(f(x),x))

报错： ‘Symbol’ object has no attribute ‘sinh’
改进：
添加一行，并把np.sinh(x)改sinh

from sympy import sinh

得到初始值
0.149

总结：方法并不好，只能拿来掌握一下odient的函数，后续使用差分法做讨论

2. 利用差分法