头歌实践教学平台 Python答案 Python 计算思维训练——常微分方程类

第1关：常微分方程求解的类实现

"""
通用一阶常微分方程：
dy/dx = f(x, y), a<=x<=b
y(x0) = y0

提供了3种数值求解方法：前向欧拉法、后向欧拉法、梯形法（改进的欧拉法）
注：后两者采用两步法进行实现
"""
 
import numpy as np

class ODE:
    def __init__(self, f, y0):  # f为微分方程函数
        self.f = f
        self.y0 = y0

    # 前向欧拉法：
    def forward(self, x):
        y = [self.y0] * len(x)
        for k in range(len(x)-1):
            h = x[k+1]-x[k]
            y[k+1] = y[k] + h * self.f(x[k], y[k])
        return y  

    # 后向欧拉法（两步）：
    def backward(self, x):
        y = [self.y0] * len(x)
           #请在下方补充代码，完成后向欧拉法
        ########### Begin ########### 
        for k in range(len(x) - 1):
            h = x[k+1] - x[k]
            y[k + 1] = y[k] + h * self.f(x[k],y[k])
            y[k + 1] = y[k] + h * self.f(x[k +1 ],y[k + 1])
        return y
        ########### End #######
        # 梯形法（两步，改进的欧拉法）：
    def trapezoidal(self, x):
        y = [self.y0] * len(x)
           #请在下方补充代码，完成梯形欧拉法
        ########### Begin ###########
        for k in range(len(x)-1):
            h = x[k + 1] - x[k]
            yp = y[k] + h * self.f(x[k],y[k])
            yc = y[k] + h * self.f(x[k] + h,yp)
            y[k + 1] = (1/2) * (yp + yc)




        
        ########### End ########### 
        return y
    
    def solve(self, x, m = 't'):
        ''' 
        ode根据方法参数m不同调用不同的方法：
        - m == 'f'或'forward'     ：前向欧拉法
        - m == 'b'或'backward'    ：后向欧拉法（两步）
        - m == 't'或'trapezoidal' ：梯形法（两步，改进的欧拉法）
        缺省：梯形法
        '''
        if m.lower() == 'f' or m.lower()=='forward':
            y = self.forward(x)
        if m.lower() == 'b' or m.lower()=='backward':
            y = self.backward(x)
        if m.lower() == 't' or m.lower()=='trapezoidal':
            y = self.trapezoidal(x)
        return np.array(y)

第2关：ODE类的继承

import numpy as np

class ODE:
    def __init__(self, f, y0):  # f为微分方程函数
        self.f = f
        self.y0 = y0

    # 前向欧拉法：
    def forward(self, x):
        y = [self.y0] * len(x)
        for k in range(len(x)-1):
            h = x[k+1]-x[k]
            y[k+1] = y[k] + h * self.f(x[k], y[k])
        return y  

    # 后向欧拉法（两步）：
    def backward(self, x):
        y = [self.y0] * len(x)
           #请在下方补充代码，完成后向欧拉法
        ########### Begin ########### 

        for k in range(len(x)-1):
            h = x[k+1]-x[k]
            y_p = y[k] + h * self.f(x[k], y[k])
            y[k+1] = y[k] + h * self.f(x[k+1], y_p)
       
        ########### End ########### 
        return y

    # 梯形法（两步，改进的欧拉法）：
    def trapezoidal(self, x):
        y = [self.y0] * len(x)
           #请在下方补充代码，完成梯形欧拉法
        ########### Begin ########### 

        for k in range(len(x)-1):
            h = x[k+1]-x[k]
            y_p = y[k] + h * self.f(x[k], y[k])
            y_c = y[k] + h * self.f(x[k+1], y_p)
            y[k+1] = 1/2 * (y_p + y_c)
        
        ########### End ########### 
        return y
    
    def solve(self, x, m = 't'):
        ''' 
        ode根据方法参数m不同调用不同的方法：
        - m == 'f'或'forward'     ：前向欧拉法
        - m == 'b'或'backward'    ：后向欧拉法（两步）
        - m == 't'或'trapezoidal' ：梯形法（两步，改进的欧拉法）
        缺省：梯形法
        '''
        if m.lower() == 'f' or m.lower()=='forward':
            y = self.forward(x)
        if m.lower() == 'b' or m.lower()=='backward':
            y = self.backward(x)
        if m.lower() == 't' or m.lower()=='trapezoidal':
            y = self.trapezoidal(x)
        return np.array(y)

class CIRCUIT(ODE):
    def f(self, t, Q):
        return (self.E - Q/self.C)/self.R

    def __init__(self, E, R, C, Q0):
           #请在下方补充代码，完成后向欧拉法
        ########### Begin ########### 
        self.E=E
        self.R=R
        self.C=C
        ODE.__init__(self,self.f,Q0)

        
        ########### End ########### 

    def __call__(self, x, m = 't'):
        return self.solve(x,m)