python 解决微分方程的操作(数值解法)

第python解决微分方程的操作(数值解法)Python求解微分方程（数值解法）

对于一些微分方程来说，数值解法对于求解具有很好的帮助，因为难以求得其原方程。

比如方程：

但是我们知道了它的初始条件，这对于我们叠代求解很有帮助，也是必须的。

那么现在我们也用Python去解决这一些问题，一般的数值解法有欧拉法、隐式梯形法等，我们也来看看这些算法对叠代的精度有什么区别？

```python

```python

importnumpyasnp

fromegrateimportodeint

frommatplotlibimportpyplotasplt

importos

#先从odeint函数直接求解微分方程

#创建欧拉法的类

classEuler:

#构造方法，当创建对象的时候，自动执行的函数

def__init__(self,h,y0):

#将对象与对象的属性绑在一起

self.h=h

self.y0=y0

self.y=y0

self.n=1/self.h

self.x=0

self.list=[1]

#欧拉法用list列表,其x用y叠加储存

self.list2=[1]

self.y1=y0

#改进欧拉法用list2列表，其x用y1叠加储存

self.list3=[1]

self.y2=y0

#隐式梯形法用list3列表，其x用y2叠加储存

#欧拉法的算法，算法返回t,x

defcountall(self):

foriinrange(int(self.n)):

y_dere=-20*self.list[i]

#欧拉法叠加量y_dere=-20*x

y_dere2=-20*self.list2[i]+0.5*400*self.h*self.list2[i]

#改进欧拉法叠加量y_dere2=-20*x(k)+0.5*400*delta_t*x(k)

y_dere3=(1-10*self.h)*self.list3[i]/(1+10*self.h)

#隐式梯形法计算y_dere3=(1-10*delta_t)*x(k)/(1+10*delta_t)

self.y+=self.h*y_dere

self.y1+=self.h*y_dere2

self.y2=y_dere3

self.list.append(float("%.10f"%self.y))

self.list2.append(float("%.10f"%self.y1))

self.list3.append(float("%.10f"%self.y2))

returnnp.linspace(0,1,int(self.n+1)),self.list,self.list2,self.list3

step=input("请输入你需要求解的步长:")

step=float(step)

work1=Euler(step,1)

ax1,ay1,ay2,ay3=work1.countall()

#画图工具plt

plt.figure(1)

plt.subplot(1,3,1)

plt.plot(ax1,ay1,'s-.',MarkerFaceColor='g')

plt.xlabel('横坐标t',fontproperties='simHei',fontsize=20)

plt.ylabel('纵坐标x',fontproperties='simHei',fontsize=20)

plt.title('欧拉法求解微分线性方程步长为'+str(step),fontproperties='simHei',fontsize=20)

plt.subplot(1,3,2)

plt.plot(ax1,ay2,'s-.',MarkerFaceColor='r')

plt.xlabel('横坐标t',fontproperties='simHei',fontsize=20)

plt.ylabel('纵坐标x',fontproperties='simHei',fontsize=20)

plt.title('改进欧拉法求解微分线性方程步长为'+str(step),fontproperties='simHei',fontsize=20)

plt.subplot(1,3,3)

plt.plot(ax1,ay3,'s-.',MarkerFaceColor='b')

plt.xlabel('横坐标t',fontproperties='simHei',fontsize=20)

plt.ylabel('纵坐标x',fontproperties='simHei',fontsize=20)

plt.title('隐式梯形法求解微分线性方程步长为'+str(step),fontproperties='simHei',fontsize=20)

plt.figure(2)

plt.plot(ax1,ay1,ax1,ay2,ax1,ay3,'s-.',MarkerSize=3)

plt.xlabel('横坐标t',fontproperties='simHei',fontsize=20)

plt.ylabel('纵坐标x',fontproperties='simHei',fontsize=20)

plt.title('三合一图像步长为'+str(step),fontproperties='simHei',fontsize=20)

ax=plt.gca()

ax.legend(('$Eular$','$fixedEular$','$trapezoid$'),loc='lowerright',title='legend')

plt.show()

os.system("pause")

对于欧拉法，它的叠代方法是：

改进欧拉法的叠代方法：

隐式梯形法：

对于不同的步长，其求解的精度也会有很大的不同，我先放一几张结果图：

补充：基于python的微分方程数值解法求解电路模型

安装环境包

安装numpy（用于调节range）和matplotlib（用于绘图）

在命令行输入

pipinstallnumpy

pipinstallmatplotlib

电路模型和微分方程

无损害，电容电压为5V，电容为0.01F，电感为0.01H的并联谐振电路

电路模型1

微分方程1

带电阻损耗的电容电压为5V，电容为0.01F，电感为0.01H的的并联谐振

电路模型2%20

微分方程2

python代码

import%20numpy%20as%20np

import%20matplotlib.pyplot%20as%20plt

L%20=%200.01%20%20#电容的值%20F

C%20=%200.01%20%20#电感的值%20L

u_0%20=%205%20%20%20#电容的初始电压

u_dot_0%20=%200

def%20equition(u,u_dot):#二阶方程

%20%20%20%20u_double_dot%20=%20-u/(L*C)

%20%20%20%20return%20u_double_dot

def%20draw_plot(time_step,time_scale):#时间步长和范围

%20%20%20%20u%20=%20u_0

%20%20%20%20u_dot%20=%20u_dot_0%20%20#初始电压和电压的一阶导数

%20%20%20%20time_list%20=%20[0]%20#时间lis

%20%20%20%20Votage%20=%20[u]%20#电压list

%20%20%20%20plt.figure()

%20%20%20%20for%20time%20in%20np.arange(0,time_scale,time_step):#使用欧拉数值计算法%20一阶近似

%20%20%20%20%20%20%20%20u_double_dot%20=%20equition(u,u_dot)%20#二阶导数

%20%20%20%20%20%20%20%20u_dot%20=%20u_dot%20+%20u_double_dot*time_step%20#一阶导数

%20%20%20%20%20%20%20%20u%20=%20u%20+%20u_dot*time_step%20#电压

%20%20%20%20%20%20%20%20time_list.append(time)%20#结果添加

%20%20%20%20%20%20%20%20Votage.append(u)%20#结果添加

%20%20%20%20%20%20%20%20print(u)

%20%20%20%20plt.plot(time_list,Votage,"b--",linewidth=1)%20#画图

%20%20%20%20plt.show()

%20%20%20%20plt.savefig("easyplot.png")

if__name__=='__main__':

draw_plot(0.0001,1)

importnumpyasnp

importmatplotlib.pyplotasplt

L=0.01#电容的值F

C=0.01#电感的值L

R=0.1#电阻值

u_0=5#电容的初始电压

u_dot_0=0

defequition(u,u_dot):#二阶方程

u_double_dot=(-R*C*u_dot-u)/(L*C)

returnu_double_dot

defdraw_plot(time_step,time_scale):#时间步长和范围

u=u_0

u_dot=u_dot_0#初始电压和电压的一阶导数

time_list=[0]#时间lis

Votage=[u]#电压list

plt.figure()

fortimeinnp.arange(0,time_scale,time_step):#使用欧拉数值计算法一阶近似

u_double_dot=equition(u,u_dot)#二阶导数

u_dot=u_dot+u_double_dot*time_step#一阶导数

u=u+u_dot*time_step#电压

time_list.append(time)#结果添加

Votage.append(u)#结果添加

print(u)

plt.plot(time_list,Votage,"b-",linewidth=1)#画图

plt.show()

plt.savefig("result