精华内容
下载资源
问答
  • 一阶常微分方程的五种解法,和二阶常微分方程的快速求解方法,方便快速回忆,根据山师老师的讲义,非常简单易推导的解法
  • 具体方程解法和背景不在赘述,见https://blog.csdn.net/Mezikov/article/details/107461970 1 代码 import math import numpy as np import matplotlib.pyplot as plt from math import e from numpy

    0 写在前面

    这篇博客是在将我上一篇博客的matlab代码移植到python中,应为后续要开展深度强化学习下的船舶减摇研究,总的来说还是在python上做这项工作比较合适。具体方程的解法和背景不在赘述,见https://blog.csdn.net/Mezikov/article/details/107461970

    1 代码

    import math
    import numpy as np
    import matplotlib.pyplot as plt
    from math import e
    from numpy import reshape
    
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    
    ##################   写在前面:这个代码只适合用于已知浪向角的规则波   ########################
    
    ts, ys, zs, ys_, zs_ = [], [], [], [], []
    # all below is coefficients of wave
    fi = math.pi / 3
    B = 1
    d = 0.341
    delta = 1.76 * 10 ** 3
    V = 1.76
    g = 9.8
    L = 6
    lamda = 1.64
    rho = 1.025 * 10 ** 3
    A_w = 5.448
    A_m = 0.346
    C_w = A_w / (L * B)
    C_p = V / (L * A_m)
    C_b = delta / (1000 * L * B * d)
    v = 4.44
    T = 0.8 * math.sqrt(lamda)
    W_0 = 2 * math.pi / T
    W_e = W_0 * (1 + v * W_0 * np.cos(fi) / g)
    x_b = -0.23
    GM_l = (L ** 2 * (5.55 * C_w + 1) ** 3 / 3450) / (C_b * d)
    M = 1.6 * 10 ** 3
    H_0 = B / (2 * d)
    a_zz = 0.8 * H_0 * C_w * M
    b_zz = (5.4 * (C_w / C_p) * (H_0 ** 0.5) - 4.7) * delta / ((g * L) ** 0.5)
    c_zz = rho * g * A_w
    a_z0 = (v / W_e ** 2) * b_zz - a_zz * x_b
    b_z0 = -b_zz * x_b + v * a_zz
    c_z0 = -rho * g * A_w * x_b
    I_00 = 0.07 * C_w * M * L ** 2
    a_00 = 0.83 * H_0 * C_p ** 2 * (0.25 * L) ** 2 * (delta / g) ** 0.5
    b_00 = 0.08 * H_0 * delta * L ** 2 / (g * L) ** 0.5
    c_00 = rho * g * (V * GM_l)
    a_0z = -(a_zz * x_b)
    b_0z = -(b_zz * x_b)
    c_0z = -rho * g * A_w * x_b - v * b_zz
    fi = math.pi / 3
    a = 0.5 * 0.17 * lamda ** (3 / 4)
    k = 2 * math.pi / lamda
    n1 = (I_00 + a_00) - (a_0z * a_z0) / (M + a_zz)
    n2 = rho * g * a * k * math.e ** (-k * d / 2) * B * d * np.sin(k * B * L * np.sin(fi) / 2) / (
                k * B * L * np.sin(fi) / 2)
    n3 = (I_00 + a_00) * (M + a_zz) - a_0z * a_z0
    
    
    #################   以上全是水动力系数的计算,下面是四阶龙格库塔算法   ######################
    def cal_F(t):
        m1 = 547 * np.cos(5.53 * t)
        m2 = 1045 * np.sin(5.53 * t)
        F_ = (I_00 + a_00) / n3 * m1 - a_z0 / n3 * m2
        return F_
    
    
    def g(y1, z1, y2, z2):
        j = (b_zz * (I_00 + a_00) - a_z0 * b_0z) / n3 * z1 + \
            (c_zz * (I_00 + a_00) - a_z0 * c_0z) / n3 * y1 + \
            (b_z0 * (I_00 + a_00) - a_z0 * b_00) / n3 * z2 + \
            (c_z0 * (I_00 + a_00) - a_z0 * c_00) / n3 * y2
        # 文献结果
        # j= 0.9716*np.cos(1.245*t-0.5566)-0.0674*np.cos(1.245*t-1.4375)-0.5252*z1-1.1424*y1-2.1741*z2-3.5594*y2
        return j
    
    
    def cal_M(t):
        m1 = 547 * np.cos(5.53 * t)
        m2 = 1045 * np.sin(5.53 * t)
        M_ = (M + a_zz) / n3 * m1 - a_0z / n3 * m2
        return M_
    
    
    def g_(y1, z1, y2, z2):
        # m1 = 547*np.cos(5.53*t)
        # m2 = 1045*np.sin(5.53*t)
        # j = (0.2580 * np.sin(5.53 * t) - 0.0042 * np.cos(5.53 * t)) - 0.2137 * z2 - 33.3415 * y2 - 0.0036 * z1 - 0.1641 * y1
        j = (b_00 * (M + a_zz) - a_0z * b_z0) / n3 * z2 + \
            (c_00 * (M + a_zz) - a_0z * c_z0) / n3 * y2 + \
            (b_0z * (M + a_zz) - a_0z * b_zz) / n3 * z1 + \
            (c_0z * (M + a_zz) - a_0z * c_zz) / n3 * y1
        # j=1.1854*np.cos(1.245*t- 1.4375)-0.006*np.cos(1.245*t-0.5566)-7.9154*z2-20.3425*y2-0.0305*z1-0.0547*y1
        return j
    
    def runge_kutta():
        """
            兄弟们,之前这里被网上的资料坑了,我发现按照oringin_try里面的办法,龙格库塔里面只有return的那个量会
            更改,其余变量都不变化,如果这是个一阶微分方程,那只有一个量在变化,可现在是二元二阶方程组啊!我现在改
            过来了,就用这种新方法了。
        """
        t = 0
        dt = 0.1
        y1 = 0
        z1 = 0
        y2 = 0
        z2 = 0
        while t <= 500:
            ys.append(y1)
            zs.append(z1)
            ys_.append(y2)
            zs_.append(z2)
            ts.append(t)
            t += dt
    
            k1 = z1  # f(t,y1,z1,y2,z2)
            k1_ = z2  # f_(t,y1,z1,y2,z2)
            l1 = cal_F(t) - g(y1, z1, y2, z2)
            l1_ = cal_M(t) - g_(y1, z1, y2, z2)
            k2 = z1 + 0.5 * dt * l1
            k2_ = z2 + 0.5 * dt * l1_
            l2 = (cal_F(t) + cal_F(t + dt)) / 2 - g(y1 + 0.5 * dt * k1, z1 + 0.5 * dt * l1, y2 + 0.5 * dt * k1_,
                                                    z2 + 0.5 * dt * l1_)
            l2_ = (cal_M(t) + cal_M(t + dt)) / 2 - g_(y1 + 0.5 * dt * k1, z1 + 0.5 * dt * l1, y2 + 0.5 * dt * k1_,
                                                      z2 + 0.5 * dt * l1_)
            k3 = z1 + 0.5 * dt * l2
            k3_ = z2 + 0.5 * dt * l2_  #  我曹之前这里z2+0.5+dt*l2,我说怎么会这么奇怪一直和Matlab对不上,写代码不能分心!
            l3 = (cal_F(t) + cal_F(t + dt)) / 2 - g(y1 + 0.5 * dt * k2, z1 + 0.5 * dt * l2, y2 + 0.5 * dt * k2_,
                                                    z2 + 0.5 * dt * l2_)
            l3_ = (cal_M(t) + cal_M(t + dt)) / 2 - g_(y1 + 0.5 * dt * k2, z1 + 0.5 * dt * l2, y2 + 0.5 * dt * k2_,
                                                      z2 + 0.5 * dt * l2_)
            k4 = z1 + dt * l3  # wo ta ma zhengde cao le 这里我为什么z1会乘以0.5啊。
            k4_ = z2 + dt * l3_  # f_(t,y1,z1,y2,z2)+0.5*dt*l3_
            l4 = cal_F(t + dt)- g(y1 + dt * k3, z1 + dt * l3, y2 + dt * k3_,
                                                    z2 + dt * l3_)  # 又乘了0.5,第四步应该都是dt而不是0.5*dt
            l4_ = cal_M(t + dt)- g_(y1 + dt * k3, z1 + dt * l3, y2 + dt * k3_, z2 + dt * l3_)
    
            y1 = y1 + dt * (k1 + 2 * k2 + 2 * k3 + k4) / 6
            z1 = z1 + dt * (l1 + 2 * l2 + 2 * l3 + l4) / 6
            y2 = y2 + dt * (k1_ + 2 * k2_ + 2 * k3_ + k4_) / 6
            z2 = z2 + dt * (l1_ + 2 * l2_ + 2 * l3_ + l4_) / 6
    
    
    ###################   主程序入口,代码到这里就算编完了  #########################3
    if __name__ == "__main__":
        runge_kutta()
        print(ys)
        # print(ys)
        plt.figure(1)
        plt.plot(ts[2000:3000], zs[2000:3000], label='埀荡')
        plt.xlabel('time(s)')
        plt.ylabel('埀荡量(m)')
        plt.legend()
        plt.show()
        plt.figure(2)
        plt.plot(ts[2000:3000], zs_[2000:3000], label='纵摇')
        plt.xlabel('time(s)')
        plt.ylabel('纵摇角(m)')
        plt.legend()
        plt.show()
        print('最大埀荡量', np.max(ys), )
        print('最大纵摇角', np.max(ys_), )
    

    2 结果

    在这里插入图片描述
    在这里插入图片描述

    3 总结

    但是不知道为什么,matlab算的和python算的结果有一点点(2%左右误差)不一样,这是为什么呢

    展开全文
  • 微分方程的数值解法主要包括两大类:有限差分法和有限单元法。这里主要介绍有限单元法。 However,对于一个只学过微积分和矩阵论的工科生来说,要了解有限元法的数学原理还是有些困难,所以这里重点是介绍有限元法...

    我们常常用微分方程来描述现实世界中的一些物理现象。由于微分方程的复杂性,只有在很简单的情况下才能得到微分方程的解析解。由于计算机的发展,采用数值方法求解微分方程的数值近似解得到广泛应用。微分方程的数值解法主要包括两大类:有限差分法和有限单元法。这里主要介绍有限单元法。

    However,对于一个只学过微积分和矩阵论的工科生来说,要了解有限元法的数学原理还是有些困难,所以这里重点是介绍有限元法的思想和具体计算方法,深层的数学原理并不涉及,主要是本人也不懂~~~

    1.微分方程基本概念

    微分方程主要分为常微分方程和偏微分方程。
    常微分方程:未知函数是单一自变量的函数,按照方程中导数的阶数,分为一阶、二阶、n阶常微分方程。
    偏微分方程:未知函数是多个自变量的函数,至少是两个自变量,同样可以分为一阶、二阶、n阶偏微分微分方程。

    二阶偏微分方程研究较多,又可以分为椭圆型方程、双曲型方程、抛物型方程等。这里不展开,具体可以查看偏微分方程数值解法相关书籍。

    2.二阶常微分方程的有限元法求解举例

    由于下面涉及大量公式,所以采用截图的方式进行介绍

    3.具体Matlab实现如下: 

    %% main program
    %一次有限元求一维常微分方程,基函数为分片线性函数
    n = 10;
    err = zeros(8, 1);
    
    % Linear interpolation
    index = zeros(8, 1);
    for i=1:8
        index(i) = 2^(i-1) * n; %10,20,40,80,160,320,640,1280
    end
    for i = 1:8
        N = index(i); %网格数量
        h = 1 / N; %网格大小
        x = 0:h:1; %节点向量
        xx = 0:1/1000:1;
        exact_solu = sin(pi.*xx); %精确解
        K = zeros(N-1); %初始化系数矩阵(刚度矩阵)
        K = K + diag(ones(N-1,1).*(2/h+2/3*h), 0) + diag((h/6-1/h).*ones(N-2,1), 1)...
              + diag((h/6-1/h).*ones(N-2,1), -1);
        F = (1+pi^2)/(h*pi^2) * (2.*sin(pi.*x(2:N))-sin(pi.*x(1:N-1))-sin(pi.*x(3:N+1)))';
        K = sparse(K);
        coeff = K\F; %Ax=b,matlab求法x=A\b
        solu = [0; coeff; 0]; %u(x)节点解
        
    %     figure()
        subplot(3,3,i)
        plot(x', solu, 'o--',xx', exact_solu, 'g--');
        legend('numerical solution', 'exact solution')
        xlabel('x')
        ylabel('U');
        xlim([-0.05,1.05])
        ylim([-0.05,1.05]);
        title(['numerical solution and exact solution N=',num2str(index(i))]);
     
    % calculate the error
        du = solu(2:N+1) - solu(1:N);
        y1 = sin(2.*pi.*x);
        dy1 = y1(2:N+1) - y1(1:N);
        y2 = sin(pi.*x);
        dy2 = y2(2:N+1) - y2(1:N);
        error = pi^2*h/2*ones(1,N) + du'.^2/h + pi/4*dy1 - 2/h*dy2.*du';
        err(i) = sqrt(sum(error));
    end
    % figure()
    subplot(3,3,9)
    plot(log2(err),'o--');
    title('误差阶')
    xlabel('N')
    ylabel('$log_{2}Error$', 'Interpreter','LaTex')
    xlim([0.8,8.2])
    
    axis on;
    set(gca,'xtick',1:1:8);
    set(gca,'xticklabel',{5,10,20,40,80,160,320,640});
    
    [order,~] = polyfit(1:1:(length(err)), log2(err)', 1);
    disp(order)
    
    

    最后计算结果如下图所示:

    参考文献:

    代码 https://blog.csdn.net/waitingwinter/article/details/106164350

    展开全文
  • 第四十九篇 二阶常微分方程的打靶法 边界值问题 当我们试图用自变量的不同值所提供的信息来解二阶或更高阶的常微分方程时,我们必须使用与前面描述的不同的数值方法。 前面提到的初值问题经常涉及到“时间”作为自...

    第四十九篇 二阶常微分方程的打靶法

    边界值问题

    当我们试图用自变量的不同值所提供的信息来解二阶或更高阶的常微分方程时,我们必须使用与前面描述的不同的数值方法。
    前面提到的初值问题经常涉及到“时间”作为自变量,解决技术要求我们按步骤“前进”,直到达到所需的解。在这种情况下,解的定义域不是有限的,因为原则上我们可以无限地沿着正方向或负方向前进。
    边值问题涉及一个有限解的区间,在这个区间内去求解。这类问题的自变量通常是空间中测量距离的坐标。一个典型的二阶边值问题可能下面的形式。
    在这里插入图片描述
    解的定义域(假设为B >a)是由x在A≤x≤B的范围内的值给出的,我们需要的是找到对应于x在这个范围内的y的值。
    本篇将考虑的数值解法分为三类
    (a)用有限差分等效代替原来的微分方程的技术。
    (b)“打靶法”,试图用一个等价初值问题代替边值问题。
    ©“加权残差”法,即猜测一个满足边界条件的试解,并调整解内的某些参数以使误差最小化。

    有限差分法

    在这种方法中,微分方程中的导数项被有限差分近似代替。正如微分方程是线性的,将会显示的那样。这个过程会导致一个线性联立方程,可以使用前面描述的方法求解。
    首先,需要定义一些有限差分近似来处理经常遇到的导数。考虑下图的解曲线,其中x轴被细分为规则的网格点,距离为h。
    在这里插入图片描述
    可以用下面任意一种方法来表示y对x的一阶导数
    在这里插入图片描述
    这些有限差分公式通过计算函数在xi附近的连接点的直线的斜率来逼近一阶导数。正向和反向公式有偏差,因为它们只包括xi的右侧或左侧的点。中心差分形式将xi左右的点考虑在内。
    高阶导数也可以用这种方法近似。例如,在xi处的二阶导数可以通过取一阶导数的逆向差来估计,如下所示
    在这里插入图片描述
    如果我们把一阶导数的前向差分公式代入,得到
    在这里插入图片描述
    类似上面可以得到yi’‘的中心差分公式
    在这里插入图片描述
    类似地,三阶导数的中心差分公式可以由
    在这里插入图片描述
    对x的四阶导数
    在这里插入图片描述
    这些高阶导数的正向和逆向差分也很容易得到,并且通过在差分公式中包含更多的点,可以获得更高的精度。
    逆向差分公式是正向差分公式的简单镜像,除了奇数导数(即,y’,y " '等),其中符号必须相反。
    一阶导数在x0处的四点正向差分公式为
    在这里插入图片描述
    逆向差分公式为
    在这里插入图片描述
    “两点边值问题的解决方案包括把x的范围分裂到n个相等的部分,每个宽度h。n = 4,如果给定的边界条件,x = A和x = B,令ξ= x0 + ih, i = 1, 2,……,其中x0 = A, xn = B。
    在每个点i = 1,2,…n−1处,微分方程以有限差分形式表示。如果微分方程是线性的,这将导致n−1个联立线性方程,在未知值yi, i = 1,2,…,n−1,其中yi表示在xi处的解。细分得越多,解的细节和精度就越高,但必须解决的联立方程也就越多。

    计算实例

    给定y’’ = 3x + 4y,边界条件y(0) = 0, y(1) = 1,利用h = 0.2有限差分求解0≤x≤1范围内的方程。
    在这里插入图片描述
    首先将微分方程写成有限差分形式。选择二阶导数公式。通常使用最低阶中心差分形式,除非有特殊的情况去使用较不准确的正向或逆向差分形式,因此
    在这里插入图片描述
    微分方程可以写成
    在这里插入图片描述
    求解域被分割成5条相等的条带,每条宽度h = 0.2。
    然后将有限差分方程写在每一个需要解的网格点上。需要时引入已知的边界条件,得到下表所示的四个方程。
    这些(非对称)方程可以用线性方程求解中的任何合适的程序求解。第二个表将数值解与精确解进行了比较
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    显然,有限差分解的精度可以通过在0≤x≤1范围内进行更多的细分来进一步提高,但代价是求解更多的方程。
    上个例子显示了一个问题,那就是如果使用了y "的高阶有限差分表示。例如,如果对y’‘使用五点中心差分公式,那么为了表示y1’‘和y4’’,就需要y在解域之外的值。为了解这个系统,还需要更多的信息,因为未知数要比方程多。
    在解域之外引入点(至少是暂时引入点)的要求经常会碰到,并且可以通过采取合并适当的边界条件来解决,如下一个示例所示。
    考虑一个边值问题
    在这里插入图片描述
    在这种情况下,y (B)是未知的,所以微分方程的有限差分需要写在x = B,即使使用最简单的中心差分公式y”,相应的“解决方案”在x5,将如下图所示。在这种情况下,我们有5个未知数,但只有4个方程。
    在这里插入图片描述
    第五个方程来自导数边界条件,它也可以写成有限差分形式,即使用中心差分
    在这里插入图片描述

    计算实例

    给定x2y’’ + 2xy’−2y = x2,边界条件y(1) = 1, y’(1.5) =−1。用h = 0.1有限差分求解1≤x≤1.5范围内的方程。
    微分方程以有限差分形式,对两个导数项使用中心差分写成如下
    在这里插入图片描述
    并应用于5个需要解的x值,如下图所示。
    得到的方程如下表所示。从下图可以看出,第五个方程在解域外引入了第六个未知数y6 = y(1.6)。中心差分形式的导数边界条件
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    用于提供第六个方程。
    这六个线性方程的解,连同精确解,给出在下表中
    在这里插入图片描述
    在这里插入图片描述
    本节给出的有限差分解在网格点相对较少的情况下表现得相当好。但并不会总是这种情况,特别是当导数在解域中发生突然的变化时。在没有精确解可供比较的情况下,建议使用两个或三个不同层次的解域来解决问题。解对网格尺寸参数h的灵敏度往往表明解的精度。
    从我们的角度来看,应该使h足够小,以达到所需的精度,但不能小于必要的值,因为这会导致过度大的方程组。

    打靶法

    打靶法试图像解决初值问题一样解决边值问题。考虑下面的二阶方程
    在这里插入图片描述
    在一个打靶法中,将解决一系列的初值问题的形式
    在这里插入图片描述
    通过有条不紊地改变初始梯度ai,我们最终可以在x = B处得到一个足够接近所需边值yB的解。有几种不同的策略来收敛所需的解,这些方法类似于寻找非线性代数方程的根。
    这里描述的方法如下图所示是选择两个初始梯度y ’ 0 (A) = a0和y ’ 1 (A) = a1,通过一步法之后,得到y0(B) = y0和y1 (B) = y₁分别在所需的边界条件yB的两边,即,y0 < yB < y1(或y1< yB < y0)。
    在这里插入图片描述
    通过线性插值,给出了初始梯度的改进估计
    在这里插入图片描述
    从而得出y(B) = y∗。
    根据下面的测试,其中一个初始梯度被∗取代,如果
    在这里插入图片描述
    然后用y代替y0和a代替a0
    或者用y代替y1和a代替a1
    在整个迭代过程中,我们的目标是保持一个高估目标边界条件的初始梯度,以及一个低估目标边界条件的初始梯度。随着迭代的进行,y∗趋于目标值,当根据准则满足收敛容差时,计算停止
    在这里插入图片描述
    迭代过程本质上就是之前描述的“假位置法”。详情可见试位法求解非线性方程

    程序如下

    程序计算的二阶常微分方程为
    在这里插入图片描述

    import numpy as np
    a0=np.zeros((2))
    k0=np.zeros((2))
    k1=np.zeros((2))
    k2=np.zeros((2))
    k3=np.zeros((2))
    y=np.zeros((2))
    nsteps=5;xa=0.0;ya=0.0;xb=1.0;yb=1.0
    a0[:]=(0.0,1.0)
    tol=0.0001;limit=25
    y0=np.zeros((nsteps+1,2))
    ystar=np.zeros((nsteps+1))
    print('--二阶常微分方程的打靶法--')
    h=(xb-xa)/nsteps
    def f74(x,y):
        f74=np.zeros((2))
        f74[0]=y[1]
        f74[1]=3.0*x**2+4.0*y[0]
        return f74
    for j in range(1,3):
        x=xa;y[0]=ya;y[1]=a0[j-1]
        for i in range(0,nsteps+1):
            y0[i,j-1]=y[0]
            k0=h*f74(x,y);k1=h*f74(x+h/2.0,y+k0/2.0)
            k2=h*f74(x+h/2.0,y+k1/2.0);k3=h*f74(x+h,y+k2)
            y=y+(k0+2.0*k1+2.0*k2+k3)/6.0;x=x+h
    if (y0[nsteps,0]-yb)*(y0[nsteps,1]-yb)>0:
        print('尝试一个新梯度')
    iters=0
    while(True):
        iters=iters+1
        astar=a0[0]+(yb-y0[nsteps,0])*(a0[1]-a0[0])/(y0[nsteps,1]-y0[nsteps,0])
        x=xa;y[0]=ya;y[1]=astar
        for i in range(0,nsteps+1):
            ystar[i]=y[0]
            k0=h*f74(x,y);k1=h*f74(x+h/2.0,y+k0/2.0)
            k2=h*f74(x+h/2.0,y+k1/2.0);k3=h*f74(x+h,y+k2)
            y=y+(k0+2.0*k1+2.0*k2+k3)/6.0;x=x+h
        if (abs(ystar[nsteps]-yb)/yb)<tol:
            print(' x              y')
            for i in range(0,nsteps+1):
                print('{:9.5e}'.format(xa+i*h),end='  ')
                print('{:9.5e}'.format(ystar[i]))
            print('迭代到收敛次数',iters)
            break
        if (ystar[nsteps]-yb)*(y0[nsteps,0]-yb)>0:
            y0[:,0]=ystar;a0[0]=astar
        else:
            y0[:,1]=ystar;a0[1]=astar
    
            
    
    

    终端输出结果如下
    在这里插入图片描述

    展开全文
  • 边值问题Matlab可用BVP4C命令,但感觉比较麻烦,下面用1stOpt求解,很简单快捷:CODE:Constant Pey=9.73, Nox=8.05, uxuy=3, bd=1, cx1e=2.3,b1=309.7,b2=2832.5, a1=27.8,a2=2.15,a3=-0.84,a4=0.935;...

    边值问题Matlab可用BVP4C命令,但感觉比较麻烦,下面用1stOpt求解,很简单快捷:CODE:

    Constant Pey=9.73, Nox=8.05, uxuy=3, bd=1, cx1e=2.3,

    b1=309.7,b2=2832.5, a1=27.8,a2=2.15,a3=-0.84,a4=0.935;

    Variable t=[0:0.025:1],x=1,y'=[0,-Pey*bd*y];

    Plot x,y,y';

    ODEFunction x'=-(Nox/b1)*(b1*x+cx1e-(a1*y+a2)/(a3*y+a4));

    y''=-Pey*bd*y'-(Nox*Pey*bd*uxuy/b2)*(b1*x+cx1e-(a1*y+a2)/(a3*y+a4));

    边值估算:

    y(t=0): 0.327786532200411

    算法: 龙格-库塔-费尔博格法(Runge-Kutta-Fehlberg Method)

    步长值: 0.025

    步长数: 40

    种群数: 5

    结果:

    t             x(t)                         y'(t)                     y(t)                         x'(t)                        y''(t)                  y'(t)

    0             1                            0                            0.327786532200411            -7.66600291732416            -24.4666416985093          0

    1             0.000693840858838046        -0.00586485658562458        0.000602760183379987        -0.00510084725252375        0.0407853057400083        -0.00586485658562458

    结果过程:

    No.        t             x(t)                         y'(t)                     y(t)                         x'(t)                        y''(t)                  y'(t)

    0        0             1                            0                            0.327786532200411            -7.66600291732416            -24.4666416985093          0

    1        0.025        0.826313919267495            -0.490075360338193           0.321182502343127            -6.27870464940025            -15.2705386363518          -0.490075360338193

    2        0.05          0.683848068843495            -0.786264397862707           0.304892497030711            -5.1579190209382             -8.81154521345653          -0.786264397862707

    3        0.075        0.566664089408245            -0.947175982341636           0.282992172618731            -4.24801699438317            -4.34185274765138          -0.947175982341636

    4        0.1           0.470051559334184            -1.01530304646892            0.258303652175927            -3.50591713326772            -1.3105072646148           -1.01530304646892

    5        0.125        0.390249601896457            -1.02128711560973            0.232742495433715            -2.89829612055369            0.68698722317898           -1.02128711560973

    6        0.15          0.324233874445947            -0.987077366676523           0.207572500259331            -2.39920298546376            1.94702073906203           -0.987077366676523

    7        0.175        0.269556313218175            -0.928258904656385           0.183592425130687            -1.98822461738467            2.68638637716003           -0.928258904656385

    8        0.2           0.224224427946598            -0.855766958939663           0.161272572460497            -1.64913802323729            3.06326085759594           -0.855766958939663

    9        0.225        0.186609425886391            -0.777149796437229           0.14085444777691             -1.36893219952477            3.19261444219091           -0.777149796437229

    10        0.25          0.155375287754656            -0.697500084306226           0.122423216362664            -1.13709376937326            3.15755303846612           -0.697500084306226

    11        0.275        0.129423296920018            -0.620141103883988           0.105960033630092            -0.94507770210612            3.01768430066738           -0.620141103883988

    12        0.3           0.10784802305656             -0.547133377656043           0.0913796756856091           -0.785907023627595           2.81532466545902           -0.547133377656043

    13        0.325        0.0899021385914759           -0.479646400522898           0.0785572032799784           -0.653865484735837           2.58009717001762           -0.479646400522898

    14        0.35          0.074968022580219            -0.418231138978627           0.0673466535604279           -0.544257391395793           2.33234938931414           -0.418231138978627

    15        0.375        0.0625347128329202           -0.363018610699521           0.057593889044556            -0.453217546785379           2.08569220018259           -0.363018610699521

    16        0.4           0.0521791123811746           -0.313863587563825           0.0491452021428162           -0.377559285865171           1.84888302113582           -0.313863587563825

    17        0.425        0.0435506113178284           -0.270447635858435           0.0418528608540996           -0.314651994050368           1.62721930221573           -0.270447635858435

    18        0.45          0.0363584698185308           -0.232352118108376           0.035578472976028            -0.262321788886875           1.42356581491271           -0.232352118108376

    19        0.475        0.0303614453595037           -0.199109100417178           0.0301948156838825           -0.218770614815521           1.23910819590991           -0.199109100417178

    20        0.5           0.0253592512680494           -0.170236100523047           0.0255866050433442           -0.182510109417003           1.07390211666811           -0.170236100523047

    21        0.525        0.0211855142304534           -0.145259103337654           0.0216505513938588           -0.152307402404205           0.927270211329725          -0.145259103337654

    22        0.55          0.0177019615432633           -0.123727136798432           0.0182949507170654           -0.127140604839182           0.798085959370831          -0.123727136798432

    23        0.575        0.0147936190504607           -0.10522084841223            0.0154389909459885           -0.106162198564426           0.684973980835008          -0.10522084841223

    24        0.6           0.0123648409053087           -0.0893568825522887          0.013011899547216            -0.0886688856491884          0.586448855258726          -0.0893568825522887

    25        0.625        0.0103360247349054           -0.0757893783829889          0.0109520199975861           -0.0740767321274094          0.501009027163578          -0.0757893783829889

    26        0.65          0.00864089210235943          -0.0642095490478649          0.00920587648845576          -0.0619006582113814          0.427198170126278          -0.0642095490478649

    27        0.675        0.00722423559559995          -0.0543440348572125          0.00772726568799934          -0.0517375017368198          0.363643217620131          -0.0543440348572125

    28        0.7           0.00604005138059095          -0.0459525242309598          0.00647639968487763          -0.0432520224132175          0.309075883307506          -0.0459525242309598

    29        0.725        0.00504999039082492          -0.0388249890925048          0.00541911382433476          -0.0361653286316663          0.262342698160537          -0.0388249890925048

    30        0.75          0.00422207307997283          -0.0327787733794859          0.00452614590083615          -0.0302453015167734          0.222407244050448          -0.0327787733794859

    31        0.775        0.00352962231729211          -0.0276556945654436          0.00377248824180515          -0.0252986667688346          0.188347254882715          -0.0276556945654436

    32        0.8           0.00295037694023563          -0.0233192611702517          0.00313681097020956          -0.0211644268962893          0.159348504512056          -0.0233192611702517

    33        0.825        0.00246575500715863          -0.0196520685790281          0.00260095269653461          -0.0177084172807521          0.134696842800273          -0.0196520685790281

    34        0.85          0.00206024116846864          -0.0165534068636294          0.00214947371190719          -0.0148187912247088          0.113769329539156          -0.0165534068636294

    35        0.875        0.00172087700369442          -0.0139370945264954          0.00176926617464773          -0.0124022733798921          0.0960251142514258        -0.0139370945264954

    36        0.9           0.00143683682371022          -0.0117295387756349          0.00144921561232049          -0.0103810490993766          0.0809964905399249        -0.0117295387756349

    37        0.925        0.00119907445004718          -0.00986801430609513        0.00117990816225777          -0.00869018040220058        0.068280395780233          -0.00986801430609513

    38        0.95          0.00100002897020513          -0.00829914727748304        0.00095337824605355          -0.00727545828585919        0.0575305147328749        -0.00829914727748304

    39        0.975        0.000833379522126888        -0.00697758825381183        0.000762891746903158        -0.00609161680467778        0.0484500671551354        -0.00697758825381183

    40        1             0.000693840858838046        -0.00586485658562458        0.000602760183379987        -0.00510084725252375        0.0407853057400083        -0.00586485658562458

    t-5927693-1-pid-8

    c1.jpg

    展开全文
  • 常微分方程解法 (四): Matlab 解法

    万次阅读 多人点赞 2019-04-30 10:33:01
    常微分方程解法求解系列博文: 常微分方程解法 (一): 常微分方程的离散化 :差商近似导数、数值积分方法、Taylor 多项式近似 常微分方程解法 (二): 欧拉(Euler)方法 常微分方程解法 (三): 龙格—库塔...
  • matlab 常微分方程数值解法 源程序代码所属分类:其他开发工具:matlab文件大小:16KB下载次数:41上传日期:2019-02-13 11:03:29上 传 者:XWLYF说明:11.1 Euler方法 38011.1.1 Euler公式的推导 38011.1.2 Euler...
  • 基于微分方程组理论和矩阵理论,采用按列比较方法和待定矩阵方法,给出了非齐次项为二次多项式与指数函数乘积的一类三维二阶常系数线性微分方程组的特解公式。对特殊情况进行了讨论,并通过算例验证了微分方程组特解...
  • 基于微分方程组理论和矩阵理论,采用待定矩阵方法和按列比较方法,给出了非齐次项为三角函数与指数函数乘积的一类三维二阶常系数线性微分方程组的特解公式,对3种特殊情况进行了讨论,并通过算例验证了微分方程组特...
  • 常微分方程数值解法——python实现

    千次阅读 2018-12-29 23:35:06
    研究生课程《应用数值分析》结课了,使用python简单记录了下常微分方程数值解法。 向前欧拉法 {yi+1=yi+hif(xi,yi)y0=y(a) \left \{ \begin{array}{lr} y_{i+1}=y_i+h_i f(x_i,y_i) \\ y_0=y(a) \end{array} \right...
  • 使用matlab解决二元二阶微分方程组的求解问题,并画出包括极坐标图在内的多幅变量间的关系图
  • 常微分方程解法 (二): 欧拉(Euler)方法

    万次阅读 多人点赞 2019-04-30 09:41:37
    上一节讲了常微分方程的三种离散化 方法:差商近似导数、数值积分、Taylor 多项式近似。 目录 §2 欧拉(Euler)方法 2.1 向前 Euler 公式、向后 Euler 公式 2.2 Euler 方法的误差估计 §3 改进的 Euler 方法 ...
  • 本节叙述常微分方程的基本理论,一阶、二阶线性常微分方程解法,以及常见的其他典型问题。
  • 常见的常微分方程的一般解法

    千次阅读 多人点赞 2020-04-30 13:41:50
    本文归纳常见的常微分方程的一般解法常微分方程,我们一般可以将其归纳为如下n类: 可分离变量的微分方程(一阶) 一阶其次(非齐次)线性微分方程(一阶) 伯努利方程(一阶) 二阶常系数微分方程(二阶) 高阶...
  • 一、线性微分方程的解的结构 1.1 二阶齐次线性方程 y′′+P(x)y′+Q(x)y=0(1) y''+P(x)y'+Q(x)y=0 \tag{1} y′′+P(x)y′+Q(x)y=0(1) 定理1:如果函数y1(x)y_1(x)y1​(x)与y2(x)y_2(x)y2​(x)是方程(1)的两个解,...
  • Matlab求常微分方程组的数值解

    万次阅读 多人点赞 2020-03-08 09:50:50
    上篇博客介绍了Matlab求解常微分方程组解析解的方法:博客地址 微分方程组复杂时,无法求出解析解时,就需要求其数值解,这里来介绍。 以下内容按照Matlab官方文档提供的方程来展开(提议多看官方文档) 介绍一下...
  • 【实例简介】和Matlab应用有关的,具体介绍常微分方程的使用和解法,原理性介绍,帮助理解。局部截断误差指的是,按()式计算由到这一步的计算值与精确值之差+。为了估计它,由展开得到的精确值是()、()两式相减(注意到...
  • matlab解常微分方程——符号解法

    千次阅读 2019-09-15 10:27:26
    1、首先得介绍一下,在matlab中解常微分方程有两种方法,一种是符号解法,另一种是数值解法。在本科阶段的微分数学题,基本上可以通过符号解法解决。 2、用matlab解决常微分问题的符号解法的关键命令是dsolve命令...
  • 2018年常微分方程MATLAB解法.doc常微分方程MATLAB解法实验四一、问题背景与实验目的求微分方程的解实际应用问题通过数学建模所归纳而得到的方程,绝大多数都是微分方程,真正能得到代数方程的机会很少.另一方面,...
  • 大三时候在跳蚤市场闲逛,从一位数学院的学长那里买了一些闲书,最近翻出来刚好有李荣华、刘播老师的《微分方程数值解法》和王仁宏老师的《数值逼近》,结合周善贵...令 ,可以将二阶微分方程转化为一阶微分方程组...
  • 第二十四讲 一阶常微分方程组

    千次阅读 2018-12-24 11:42:07
    一,常微分方程组: 必须满足:自变量只有1个,因变量有多个 一阶形式:,,x和y是因变量,t是自变量 二,什么样的常微分方程组是线性方程组? 必须满足:x和y以线性形式出现 例如:,,a、b、c和d可以是t的函数...
  • 常微分方程的数值解法及仿真

    热门讨论 2010-01-07 23:40:30
    三、 一阶常微分方程组的数值解法 2 四、 仿真算例 4 仿真1 应用欧拉法 4 仿真2 应用二阶龙格-库塔法 5 仿真3 应用四阶龙格-库塔法 6 附录 Matlab程序 7 1. 欧拉法程序 7 2. 二阶龙格-库塔法程序 8 3. 四阶龙格-库塔...
  • 高阶线性微分方程-常微分方程

    千次阅读 2017-09-01 15:15:00
    高阶线性微分方程-常微分方程 这里讨论常微分方程常微分方程的阶数就是函数求导的最高次数。这里以二阶线性微分方程为例。 形如方程5的称为二阶线性微分方程。   线性的概念定义为:   下面讨论 二阶...
  • 常微分方程的数值解法 一阶常微分方程数值解的C语言编程实现导读:就爱阅读网友为您分享以下“一阶常微分方程数值解的C语言编程实现”资讯,希望对您有所帮助,感谢您对92的支持!一阶常微分方程数值解的C语言编程...
  • 本文主要介绍matlab中求解常微分方程(组)的dsolve和ode系列...dsolve函数用于求常微分方程组的精确解,也称为常微分方程的符号解。如果没有初始条件或边界条件,则求出通解;如果有,则求出特解。 1)函数格式  ...
  • 本文主要介绍matlab中求解常微分方程(组)的dsolve和ode系列函数,并通过例子加深读者的理解。 一、符号介绍    D: 微分符号;D2表示二阶微分,D3表示...dsolve函数用于求常微分方程组的精确解,也称为常...
  • 使用龙格库塔法求解二阶微分方程。可以设置仿真步长、初值,轻松更改函数
  • Matlab求常微分方程组的解析解

    千次阅读 2020-03-07 21:53:14
    Matlab求常微分方程的解析解 最近同学毕设需要求解循坏摆的微分方程,我在帮忙过程中学习了一下常微分方程的解析解和数值解的求法,在此分享。 以下讲解遵循Matlab官方文档提供的方程和写法。 一阶常微分方程求解...
  • 所谓数值解法,就是设法将常微分方程离散化,建立差分方程,给出解在一些离散点上的近似值.问题 7.1 一阶常微分方程初值问题的一般形式{y′=f(x,y),a⩽x⩽by(a)=α\begin{equation} \left \{ \begin{aligned} & y'=f(x,...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 659
精华内容 263
热门标签
关键字:

二阶常微分方程组解法