精华内容
下载资源
问答
  • 用Eular法解常微分方程组数值解,使用了细胞数组,代码简洁,除注释外的有效代码只有二十行左右。(几年前上传的程序了,当时要20积分,现在为大家降到5个积分)
  • 基于python,采用四阶龙格-库塔算法进行常微分方程组数值求解,最终用Apollo卫星运动轨迹作为示例

    预备知识

    包括最常用的四阶Ronge-Kutta数值方法以及四阶Adams预估-校正格式

    四阶R-K

    之所以是四阶R-K,是因为三阶精度太低,在步长略大时无法满足正常求解精度要求,而五阶以上虽然精度很高,但算法耗时。四阶R-K较为折中,既可满足精度要求,对于复杂计算其计算速度也可接受。
    四阶R-K公式的标准格式如下:
    对于形如   d y d x = f ( x , y ) \ \frac{dy}{dx}=f(x,y)  dxdy=f(x,y) 的常微分方程:

      y i + 1 = y i + h 6 ( K 1 + 2 K 2 + 2 K 3 + K 4 )   K 1 = f ( x i , y i )   K 2 = f ( x i + h 2 , y i + h 2 K 1 )   K 3 = f ( x i + h 2 , y i + h 2 K 2 )   K 4 = f ( x i + h , y i + h K 3 ) (1) \ y_{i+1}=y_i+\frac{h}{6}(K_1+2K_2+2K_3+K_4)\\ \ K_1=f(x_i,y_i)\\ \ K_2=f(x_i+\frac{h}{2},y_i+\frac{h}{2}K_1)\\ \tag{1} \ K_3=f(x_i+\frac{h}{2},y_i+\frac{h}{2}K_2)\\ \ K_4=f(x_i+h,y_i+hK_3)  yi+1=yi+6h(K1+2K2+2K3+K4) K1=f(xi,yi) K2=f(xi+2h,yi+2hK1) K3=f(xi+2h,yi+2hK2) K4=f(xi+h,yi+hK3)(1)

    对于常微分方程组来说,将   y \ y  y   f \ f  f 向量化即可,其中   h \ h  h 为步长,应尽量取小。

    四阶Adams预估-校正公式

    四阶Adams预估-校正格式可以看做一种改进的Euler格式,其计算效率相较四阶R-K更高,但其精度不如四阶R-K,在各种工程应用中应做出各种尝试后选取适合实际的算法。
    四阶Adams预估-校正公式标准格式如下:
    对于形如   d y d x = f ( x , y ) \ \frac{dy}{dx}=f(x,y)  dxdy=f(x,y) 的常微分方程:

    预估:   y ˉ i + 1 = y i + h 24 ( 55 f i − 59 f i − 1 + 37 f i − 2 − 9 f i − 3 ) \ \bar{y}_{i+1}=y_i+\frac{h}{24}(55f_i-59f_{i-1}+37f_{i-2}-9f_{i-3})  yˉi+1=yi+24h(55fi59fi1+37fi29fi3)
    校正:   y i + 1 = h 24 ( 9 f ( x i + 1 , y ˉ i + 1 ) + 19 f i − 5 f i − 1 + f i − 2 ) (2) \ y_{i+1}=\frac{h}{24}(9f(x_{i+1},\bar{y}_{i+1})+19f_i-5f_{i-1}+f_{i-2}) \tag{2}  yi+1=24h(9f(xi+1,yˉi+1)+19fi5fi1+fi2)(2)
    其中   i = 3 , 4 , 5... , n − 1 \ i=3,4,5...,n-1  i=3,4,5...,n1
    由于在预估式中出现了   f i − 3 , f i − 2 , f i − 1 \ f_{i-3},f_{i-2},f_{i-1}  fi3,fi2,fi1,因此在启动Adams预估校正算法之前需要用其他算法进行前三个时间步的运算,较为普遍的是采用四阶R-K公式进行启动计算。

    实战演练

    理论推导

    本文采用难度中等的常微分方程组进行演练,其物理含义为Apollo卫星的运动轨迹   ( x , y ) \ (x,y)  (x,y)
    其满足以下方程组:
      x ¨ = 2 y ˙ + x − μ ∗ ( x + μ ) r 1 3 − μ ( x + μ ∗ ) r 2 3 y ¨ = − 2 x ˙ + y − μ ∗ y r 1 3 − μ y r 2 3 (2) \ \ddot{x}=2\dot{y}+x-\frac{\mu^*(x+\mu)}{{r_1}^3}-\frac{\mu(x+\mu^*)}{{r_2}^3}\\ \tag{2} \ddot{y}=-2\dot{x}+y-\frac{\mu^*y}{{r_1}^3}-\frac{\mu y}{{r_2}^3}  x¨=2y˙+xr13μ(x+μ)r23μ(x+μ)y¨=2x˙+yr13μyr23μy(2)
    其中   μ = 1 / 82.45 , μ ∗ = 1 − μ r 1 = ( x + μ ) 2 + y 2 , r 2 = ( x − μ ∗ ) 2 + y 2 \ \mu=1/82.45, \mu^*=1-\mu\\ r_1=\sqrt{(x+\mu)^2+y^2},r_2=\sqrt{(x-\mu^*)^2+y^2}  μ=1/82.45,μ=1μr1=(x+μ)2+y2 ,r2=(xμ)2+y2
    之后对常微分方程组进行变量代换,取   x 1 = x , x = x ˙ , x 3 = y , x 4 = y ˙ \ x_1=x,x=\dot{x},x_3=y,x_4=\dot{y}  x1=x,x=x˙,x3=y,x4=y˙
    则式   ( 2 ) \ (2)  (2)可化为下式   ( 3 ) \ (3)  (3)
      x 1 ˙ = x 2   x 2 ˙ = 2 x 4 + x 1 − μ ∗ ( x 1 + μ ) / r 1 3 − μ ( x 1 − μ ∗ ) / r 2 3   x 3 ˙ = x 4   x 4 ˙ = − 2 x 2 + x 3 − μ ∗ x 3 / r 1 3 − μ x 3 / r 2 3 (3) \ \dot{x_1}=x_2\\ \ \dot{x_2}=2x_4+x_1-\mu^*(x_1+\mu)/{r_1}^3-\mu(x_1-\mu^*)/{r_2}^3\\ \ \dot{x_3}=x_4\\ \ \dot{x_4}=-2x_2+x_3-\mu^*x_3/{r_1}^3-\mu x_3/{r_2}^3 \tag{3}  x1˙=x2 x2˙=2x4+x1μ(x1+μ)/r13μ(x1μ)/r23 x3˙=x4 x4˙=2x2+x3μx3/r13μx3/r23(3)
    同时由于式   ( 2 ) \ (2)  (2)的初始条件为
      x ( 0 ) = 1.2 , y ( 0 ) = 0 , x ˙ ( 0 ) = 0 , y ˙ ( 0 ) = − 1.04935751 \ x(0)=1.2,y(0)=0,\dot{x}(0)=0,\dot{y}(0)=-1.04935751  x(0)=1.2,y(0)=0,x˙(0)=0,y˙(0)=1.04935751
    此时我们得到了常微分方程组,初始条件,利用四阶R-K公式,设定合适的步长即可对其进行求解。

    python实现

    下面对其python代码实现进行展示:

    创建求解常微分方程组的简单类

    # -*- coding: utf-8 -*-
    """
    Created on Mon Jan 25 11:12:42 2021
    
    @author: xinghy
    """
    # 1.2  0  0  -1.04935751'''initial'''
    import numpy as np
    
    class Ode_crew():
        '''equation crew'''
        def __init__(self, x1, x2, x3, x4, t_initial, t_last, t_iter):    
            self.x1 = [x1]
            self.x2 = [x2]
            self.x3 = [x3]
            self.x4 = [x4]
            self.time_pool = [t_initial]
            self.t_latest = t_initial
            self.t_iter = t_iter
            self.t_last = t_last
            self.x = [self.x1, self.x2, self.x3, self.x4]
            
        def f(self, x1, x2, x3, x4):
        '''方程组'''
            x1_ = x2
            x2_ = 2 * x4 + x1 - (1 - 1/82.45)*(x1 + 1/82.45)/np.sqrt((x1 + 1/82.45)**2+x3**2)**3 - (
                1/82.45)*(x1 -1 + 1/82.45)/np.sqrt((x1 - 1 + 1/82.45)**2 + x3**2)**3
            x3_ = x4
            x4_ = -2 * x2 + x3 - (1 - 1/82.45)*x3/np.sqrt((x1 + 1/82.45)**2+x3**2)**3 - (
                1/82.45)*x3/np.sqrt((x1- 1 + 1/82.45)**2 + x3**2)**3
            f = [x1_, x2_, x3_, x4_]
            return f
        
        def four_RK(self):
        '''四阶R-K'''
            K1 = self.f(self.x1[-1], self.x2[-1], self.x3[-1], self.x4[-1])
            K2 = self.f(self.x1[-1] + self.t_iter/2*K1[0], self.x2[-1] + self.t_iter/2*K1[1],
                        self.x3[-1] + self.t_iter/2*K1[2], self.x4[-1] + self.t_iter/2*K1[3])
            K3 = self.f(self.x1[-1] + self.t_iter/2*K2[0], self.x2[-1] + self.t_iter/2*K2[1],
                        self.x3[-1] + self.t_iter/2*K2[2], self.x4[-1] + self.t_iter/2*K2[3])
            K4 = self.f(self.x1[-1] + self.t_iter*K3[0], self.x2[-1] + self.t_iter*K3[1],
                        self.x3[-1] + self.t_iter*K3[2], self.x4[-1] + self.t_iter*K3[3])
            x_ori = np.array([self.x1[-1], self.x2[-1], self.x3[-1], self.x4[-1]])
            x_new = x_ori + np.array(self.t_iter/6 * (np.array(K1) + 2*np.array(K2) + 2*
                                                      np.array(K3) + np.array(K4)))
            self.x1.append(x_new[0])
            self.x2.append(x_new[1])
            self.x3.append(x_new[2])
            self.x4.append(x_new[3])
            self.t_latest += self.t_iter
            self.time_pool.append(self.t_latest)
    

    之后将各种条件代入即可用指定算法进行运算:

    import numpy as np
    import matplotlib.pyplot as plt
    
    from ode import Ode_crew
    
    odecrew = Ode_crew(1.2, 0, 0, -1.04935751, 0, 20, 0.001)
    while odecrew.t_latest < odecrew.t_last - 0.0001*0.1:
        odecrew.four_RK()
        # print(odecrew.time_pool[-1])
    x_ = odecrew.x[0]
    y_ = odecrew.x[2]
    plt.scatter(x_, y_, c=y_, cmap=plt.cm.Blues, s=0.2)
    plt.title('Apollo trajectory', fontsize=25)
    plt.xlabel('x', fontsize=10)
    plt.ylabel('y', fontsize=10)
    

    运行后即可得到如下图所示的Apollo运行轨迹图:
    由图可以看到卫星的运行轨迹结果,同时添加了一些颜色渐变的效果,比较简陋,之后有时间可以精进一下。

    附录

    前面也提到了四阶Adams预估-校正公式,是一个与四阶R-K算法总是同时出现的常用算法,但笔者运算之后发现四阶Adams所需步长要远小于四阶R-K所需步长,其精度总是与四阶R-K差几个量级,可能是方程组的固有性质所致,也可能是四阶Adams的固有弊端,在这里就不展开Adams的结果了~

    以上,请各位巨巨批评指正!

    展开全文
  • MATLAB常微分方程组数值求解方程与方程组的数值解PPT学习教案.pptx
  • 主函数比较% 根据文献自己编写q = 0.5; %分数阶阶数fdefun = @(t,y) [2/gamma(3-q)*t^... %一元微分方程y0 =0;%初值y0为列向量h = 2^(-6);%步长tspan = [0,2];[t,y] = fdewfc(q,fdefun,tspan,y0,h);figure(1)plot(t...

    f65315825342717bc64609575d86d66e.png

    主函数比较

    % 根据文献自己编写

    q = 0.5; %分数阶阶数

    fdefun = @(t,y) [2/gamma(3-q)*t^(2-q)-1/gamma(2-q)*t^(1-q)-y+t^2-t]; %一元微分方程

    y0 =0;%初值y0为列向量

    h = 2^(-6);%步长

    tspan = [0,2];

    [t,y] = fdewfc(q,fdefun,tspan,y0,h);

    figure(1)

    plot(t,y(1,1:end)) ;

    xlabel('t'); ylabel('y(t)');

    hold on

    plot(t,t.^2-t+0.5,'r-.')%解析解

    %plot(t,-2./(t+1),'r')

    %%%%%%与意大利 Roberto Garrappa, University of Bari, Italy 结果比较

    t0=tspan(1); tfinal = tspan(2);

    [t, y_fde12] = fde12(q,fdefun,t0,tfinal,y0,h);

    plot(t,y_fde12+1,'k.')

    Jfdefun = @(t,y)[-1];%Jacobi阵 此处退化为导数

    [t, y_flmm2] = flmm2(q,fdefun,Jfdefun,t0,tfinal,y0,h);

    plot(t,y_flmm2+1.5,'m--')

    legend('本方法','解析解','Roberto Garrappa方法1','Roberto Garrappa方法2','location','northwest') ;

    title('预估-校正ABM法与其他方法比较');

    子函数

    function [t,y] = fdewfc(q,fdefun,tspan,y0,h)

    % 根据文献自己编写

    %q = 0.5;

    %fdefun = @(t,y) [2/gamma(3-q)*t^(2-q)-1/gamma(2-q)*t^(1-q)-y+t^2-t]; %一元微分方程

    %初值y0为列向量

    %步长h = 2^(-6) ;

    t0 = tspan(1) ; tfinal = tspan(2);

    t = t0:h:tfinal;% 时间向量

    itn = length(t)-1;%ceil((tfinal-t0)/h);

    rd = length(y0); %向量维数

    y = zeros(rd,itn+1); y(:,1) = y0; %赋初值 按列向量排

    gq1 = 1./gamma(q); %预估用常数

    gq2 = (h.^q)./gamma(q+2);%校正用常数

    %%%%%%%%%%主循环

    for n = 0:itn-1

    % 生成B b_{j,n+1} = h^{\alpha}/\alpha((n+1-j)^{\alpha}-(n-j)^{\alpha})

    % 共n+1项,这里用 q 表示 \alpha

    jv=0:n; B = h^q/q*((n+1-jv).^q-(n-jv).^q);% n+1个元素;必须注意:每步迭代更新B

    % 生成A a_{0,n+1} = n^{\alpha+1}-(n-\alpha)(n+1)^{\alpha}, j=0,

    %       a_{j,n+1} =

    %       (n-j+2)^{\alpha+1}+(n-j)^{\alpha+1}-2(n-j+1)^{\alpha+1},1\leqslant

    %       j\leqslant n,

    %       a_{n+1,n+1}=1,j=n+1.

    % 共n+1项,这里用 q 表示 \alpha

    jv = 1:n; A = [n^(q+1)-(n-q)*(n+1)^q,(n-jv+2).^(q+1)+(n-jv).^(q+1)-2*(n-jv+1).^(q+1)]; % n+1个元素;每步迭代更新 A

    %%%% 预估 yp

    fdefunp = y0; %B = fdeB(q,h,n);%必须注意:每步迭代更新B

    for j = 1:n+1

    fdefunp = fdefunp + gq1.*B(:,j).*fdefun(t(j),y(:,j));

    end

    %%%%后一个子循环求和 %第1列和第3列合并

    fdefunsumn = 0; %fdefunsumn 应为与y0同型的零列向量,但此处写为0,无妨

    %A = fdeA(q,n); %必须注意:每步迭代更新A

    for j = 1:n+1

    fdefunsumn = fdefunsumn + A(:,j).*fdefun(t(j),y(:,j));%0~~n项

    end

    %%%% 校正

    y(:,n+2) = y0 + gq2.*(fdefun(t(n+2),fdefunp)+ fdefunsumn); %n

    end

    转载本文请联系原作者获取授权,同时请注明本文来自王福昌科学网博客。

    链接地址:http://blog.sciencenet.cn/blog-292361-1079208.html

    上一篇:用MATLAB模拟简单小球跳动

    下一篇:怀柔喇叭沟之秋随拍

    展开全文
  • 常微分方程组数值求解;主要内容 数值求解常微分方程组函数概述 非刚性/刚性常微分方程问题求解 隐式微分方程组求解 微分代数方程(DAE)与延迟微分方程(DDE)求解 边值问题求解 ;第一节数值求解常微分方程组函数概述;一...
  • Matlab求常微分方程组数值解

    万次阅读 多人点赞 2020-03-08 09:50:50
    上篇博客介绍了Matlab求解常微分方程组解析解的方法:博客地址 微分方程组复杂时,无法求出解析解时,就需要求其数值解,这里来介绍。 以下内容按照Matlab官方文档提供的方程来展开(提议多看官方文档) 介绍一下...

    上篇博客介绍了Matlab求解常微分方程组解析解的方法:博客地址
    微分方程组复杂时,无法求出解析解时,就需要求其数值解,这里来介绍。
    以下内容按照Matlab官方文档提供的方程来展开(提议多看官方文档)

    介绍一下核心函数ode45()

    一般形式:[t,y] = ode45(odefun,tspan,y0) 其中 tspan = [t0 tf]
    功能介绍:求微分方程组 y′=f(t,y) 从 t0 到 tf 的积分,初始条件为 y0。解数组 y 中的每一行都与列向量 t 中返回的值相对应。

    1. 一阶微分方程求解(简单调用即可)

    方程:y’=2*t
    代码

    tspan=[1 6]; %定义自变量x的取值空间为1-6
    y0=0;%定义因变量的初值,当x=1(x取值空间的第一个数)时,y0=0
    [t,y]=ode45(@(t,y) 2*t,tspan,y0); %定义函数y'=2*t,使用ode45求解
    plot(t,y,'-o'); %绘制求得的数值曲线
    

    说明:简单的odefun参数就是这个形式,@(x,y) fun 中间是空格,不能用逗号
    结果
    在这里插入图片描述

    2. 二阶微分方程求解(引入函数文件)

    方程:范德波尔方程 y1’’-u(1-y1²)*y1’+y1=0;(这里设u=1)
    代码
    定义输入的方程,以函数形式定义

    function dydt=odefun(t,y)
    %二阶方程为y1''-(1-y1²)*y1'+y1=0;
    %降阶为两个方程:y1'=y2; 
    %               y2'=(1-y1²)*y2-y1;
    %t虽然没有使用,但必须要作为参数写入
    dydt=[y(2);(1-y(1)^2)*y(2)-y(1)];
    

    求解作图

    tspan=[0 20]; %定义自变量x的取值空间为0-20
    y0=[2;0];%定义因变量的初值,当x=0时,y1=2,y2=y1'=0;
    [t,y]=ode45(@odefun,tspan,y0); %使用ode45求解
    %%下面为作图过程,不解释
    plot(t,y(:,1),'-o',t,y(:,2),'-o')
    title('Solution of van der Pol Equation (\mu = 1) with ODE45');
    xlabel('Time t');
    ylabel('Solution y');
    legend('y_1','y_2')
    

    说明:先进行降阶,使y2=y1’,这时候y2’=y1’’,然后把两个方程写成列向量形式就行了
    结果
    在这里插入图片描述

    3. 求解微分方程组(和2类似)

    这里就和求解二阶方程类似的,只不过不需要降阶,仍旧需要一个函数来定义方程组。我们这里不用官方文档的例子,用同学的循坏摆问题来进行演示。
    方程在这里插入图片描述
    给定的初值(w接近0,但实际上不能设置为0):
    在这里插入图片描述
    代码
    定义输入的方程

    function dRvw=func(t,Rvw)
    %% 函数功能:为ode45提供微分方程
    %输入:t:时间序列,就是θ;Rvw:因变量,Rvw(1)代表R,Rvw(2)代表v,Rvw(3)代表w
    %输出:dRvw:因变量的一阶微分,dRvw(1)代表dR,dRvw(2)代表dv,dRvw(3)代表dw
    %% 初始化因变量的一阶微分,3×1的向量
    dRvw=zeros(3,1);
    %% 参数初始化
    r=0.01;u=0.1;g=9.8;M=10;m=1;
    %% 输入微分方程式
    dRvw(1)=-Rvw(2)/Rvw(3)-r;
    dRvw(2)=(M*g-(m*Rvw(3)^2*Rvw(1)+m*g*sin(t)) *exp(u*(t+pi/2)))/(Rvw(3)*(M+m));
    dRvw(3)=Rvw(3)*r/Rvw(1)+2*Rvw(2)/Rvw(1)+g*cos(t)/(Rvw(3)*Rvw(1));
    

    求解作图

    %% 参数初始化 start_Theta是θ的初始值 end_Theta是θ的结束值
    %R是半径初值;v是线速度初值;w是角速度初值
    start_Theta=0;end_Theta=2*pi;R=1;v=0;w=1e-5;
    %% 使用ode45方法计算微分方程组func的数值解
    %func是带有方程组的函数 
    %[start_Theta end_Theta]是自变量范围
    %[R;v;w]是方程初值
    %T是自变量的数组,Rvw是对应的因变量的数值。Rvw(:,1)=R;Rvw(:,2)=v;Rvw(:,3)=w;
    [T,Rvw]=ode45(@func,[start_Theta end_Theta],[R;v;w]);
    %% 组合自变量和因变量。TRvw(:,1)=θ;TRvw(:,2)=R;TRvw(:,3)=v;Trvw(:,4);
    TRvw=[T,Rvw];
    %% 绘制自变量-因变量图像.figure1是R-θ,figure2是v-θ,figure3是w-θ
    plot(T,Rvw(:,1),'-o',T,Rvw(:,2),'-o',T,Rvw(:,3),'-o')
    title('自变量-因变量图像');
    xlabel('Angle θ');
    ylabel('Solution');
    legend('R','v','w')
    

    说明:注释的应该是比较清楚的,把三个方程写成列向量的形式就行

    PS:有些人和我说不能运行,然后我看了他们出错的原因,有点儿哭笑不得。出错的基本上都是运行上面的dRvw=func(t,Rvw)这个函数的。说明一下,这是有参数的函数,不给参数不能直接运行的。下面的求解作图脚本才是需要运行的哈,它调用了函数,才得到的结果。如果大家还发现什么问题,欢迎私戳或评论。

    PS+ 有了PS之后,还是很多人问我参数的问题,我在这里直接把文件给大家:cupt.zip
    提取码:6k8n。
    解压后,有两个文件,大家直接运行cupt.m即可。

    结果
    在这里插入图片描述

    4. 更多形式

    讲到这里,大部分我们用到的微分方程形式都可以求解了,Matlab还支持带有时变项和额外参数的微分方程求解,这里不再赘述,大家可以自行参阅官方文档。

    展开全文
  • 主函数比较% 根据文献自己编写q = 0.5; %分数阶阶数fdefun = @(t,y) [2/gamma(3-q)*t^... %一元微分方程y0 =0;%初值y0为列向量h = 2^(-6);%步长tspan = [0,2];[t,y] = fdewfc(q,fdefun,tspan,y0,h);figure(1)plot(t...

    f65315825342717bc64609575d86d66e.png

    主函数比较

    % 根据文献自己编写

    q = 0.5; %分数阶阶数

    fdefun = @(t,y) [2/gamma(3-q)*t^(2-q)-1/gamma(2-q)*t^(1-q)-y+t^2-t]; %一元微分方程

    y0 =0;%初值y0为列向量

    h = 2^(-6);%步长

    tspan = [0,2];

    [t,y] = fdewfc(q,fdefun,tspan,y0,h);

    figure(1)

    plot(t,y(1,1:end)) ;

    xlabel('t'); ylabel('y(t)');

    hold on

    plot(t,t.^2-t+0.5,'r-.')%解析解

    %plot(t,-2./(t+1),'r')

    %%%%%%与意大利 Roberto Garrappa, University of Bari, Italy 结果比较

    t0=tspan(1); tfinal = tspan(2);

    [t, y_fde12] = fde12(q,fdefun,t0,tfinal,y0,h);

    plot(t,y_fde12+1,'k.')

    Jfdefun = @(t,y)[-1];%Jacobi阵 此处退化为导数

    [t, y_flmm2] = flmm2(q,fdefun,Jfdefun,t0,tfinal,y0,h);

    plot(t,y_flmm2+1.5,'m--')

    legend('本方法','解析解','Roberto Garrappa方法1','Roberto Garrappa方法2','location','northwest') ;

    title('预估-校正ABM法与其他方法比较');

    子函数

    function [t,y] = fdewfc(q,fdefun,tspan,y0,h)

    % 根据文献自己编写

    %q = 0.5;

    %fdefun = @(t,y) [2/gamma(3-q)*t^(2-q)-1/gamma(2-q)*t^(1-q)-y+t^2-t]; %一元微分方程

    %初值y0为列向量

    %步长h = 2^(-6) ;

    t0 = tspan(1) ; tfinal = tspan(2);

    t = t0:h:tfinal;% 时间向量

    itn = length(t)-1;%ceil((tfinal-t0)/h);

    rd = length(y0); %向量维数

    y = zeros(rd,itn+1); y(:,1) = y0; %赋初值 按列向量排

    gq1 = 1./gamma(q); %预估用常数

    gq2 = (h.^q)./gamma(q+2);%校正用常数

    %%%%%%%%%%主循环

    for n = 0:itn-1

    % 生成B b_{j,n+1} = h^{\alpha}/\alpha((n+1-j)^{\alpha}-(n-j)^{\alpha})

    % 共n+1项,这里用 q 表示 \alpha

    jv=0:n; B = h^q/q*((n+1-jv).^q-(n-jv).^q);% n+1个元素;必须注意:每步迭代更新B

    % 生成A a_{0,n+1} = n^{\alpha+1}-(n-\alpha)(n+1)^{\alpha}, j=0,

    %       a_{j,n+1} =

    %       (n-j+2)^{\alpha+1}+(n-j)^{\alpha+1}-2(n-j+1)^{\alpha+1},1\leqslant

    %       j\leqslant n,

    %       a_{n+1,n+1}=1,j=n+1.

    % 共n+1项,这里用 q 表示 \alpha

    jv = 1:n; A = [n^(q+1)-(n-q)*(n+1)^q,(n-jv+2).^(q+1)+(n-jv).^(q+1)-2*(n-jv+1).^(q+1)]; % n+1个元素;每步迭代更新 A

    %%%% 预估 yp

    fdefunp = y0; %B = fdeB(q,h,n);%必须注意:每步迭代更新B

    for j = 1:n+1

    fdefunp = fdefunp + gq1.*B(:,j).*fdefun(t(j),y(:,j));

    end

    %%%%后一个子循环求和 %第1列和第3列合并

    fdefunsumn = 0; %fdefunsumn 应为与y0同型的零列向量,但此处写为0,无妨

    %A = fdeA(q,n); %必须注意:每步迭代更新A

    for j = 1:n+1

    fdefunsumn = fdefunsumn + A(:,j).*fdefun(t(j),y(:,j));%0~~n项

    end

    %%%% 校正

    y(:,n+2) = y0 + gq2.*(fdefun(t(n+2),fdefunp)+ fdefunsumn); %n

    end

    转载本文请联系原作者获取授权,同时请注明本文来自王福昌科学网博客。

    链接地址:http://blog.sciencenet.cn/blog-292361-1079208.html

    上一篇:用MATLAB模拟简单小球跳动

    下一篇:怀柔喇叭沟之秋随拍

    展开全文
  • 常微分方程的数值解法 一阶常微分方程数值解的C语言编程实现导读:就爱阅读网友为您分享以下“一阶常微分方程数值解的C语言编程实现”资讯,希望对您有所帮助,感谢您对92的支持!一阶常微分方程数值解的C语言编程...
  • 用改进欧拉法求一阶常微分方程数值解,计算结果精确
  • 求解常微分方程组 依赖 用 Fortran 90 编写的代码 gfortran 编译器 使用 Matlab/Octave 绘制解决方案 如何使用 运行代码 代码在 Fortran 90 中运行,您将需要一个 Fortran 编译器,例如 gfortran。 在代码中更改了...
  • 常微分方程数值解,MATLAB的命令格式为: [t,y]=solver('odefun',tspan,y0,options)其中solver选择ode45等函数名,odefun为根据待解方程或方程组编写的m文件名,tspan为自变量的区间[t0,tf],即准备在那个...
  • 求解一阶常微分方程数值方法(单步和多步)。 方法包括: 1.欧拉方法2. 亨氏法3. 四阶 Runge Kutta 方法4. Adams-Bashforth 方法5. Adams-Moulton 方法这些方法通常用于求解 IVP,一阶初始值问题 (IVP) 被定义为...
  • *题目描述:*求常微分方程y’=y-2x/y &...常微分方程数值解 1、y'=y-2x/y && y(0)=1; 2、解析解为y=sqrt(1+2x) */ #include <stdio.h> #include <math.h> #define N 5 //离散点个数,分别取5,
  • matlab 常微分方程数值解法 源程序代码所属分类:其他开发工具:matlab文件大小:16KB下载次数:41上传日期:2019-02-13 11:03:29上 传 者:XWLYF说明:11.1 Euler方法 38011.1.1 Euler公式的推导 38011.1.2 Euler...
  • manlab软件应用试验题目专业 序号 姓名 日期实验3 常微分方程数值解【实验目的】1.掌握用MATLAB求微分方程初值问题数值解的方法;2.通过实例学习微分方程模型解决简化的实际问题;3.了解欧拉方法和龙格库塔方法的...
  • 本文是自己写的关于怎样利用MATLAB求解常微分方程数值解的,文中从Euler法讲起,最后总结了常用的odeXX的用法及其原理,其中包含各个函数怎样使用的MATLAB代码
  • 用ODE(ordinary differential equations)23和ODE45求解微分方程 ODE(ordinary differential equations)23 通过使用二阶和三阶隆哥-库塔方法积分微分方程得出
  • 求解的算法,编写成熟上机调试出结果,要求所编程序适用于任何一阶常微分方程数值解问题,即能解决这一类问题,而不是某一个问题。 试验中以下列数据验证程序的正确性。 求 y'=-xy^2 y(0)=2 (0) 步长h=0.25 效果:...
  • 数值分析课程设计常微分方程组初值问题数值解的实现和算法分析毕业论文.doc
  • 常微分函数表达式: import scipy.integrate as spi import numpy as np from matplotlib import pyplot as plt import warnings warnings.filterwarnings(‘ignore’) 参数 ATR=0.2 ATAR=0.3 SR=0.1 SAR=0.3 SRC=...
  • 我是SymPy和Python的新手,我目前正在使用Python 2.7和SymPy 0.7.5,其目标是:a)从文本文件中读取微分方程组b)解决系统问题我已经读过this question和this other question,它们几乎就是我要找的,但我还有一个问题...
  • 常微分数值解matlab代码数值分析 两个分区: MATLAB 一些线性代数代码 C++ 求解非线性方程 求解方程组 插值和曲线拟合 函数逼近 数值微分和积分 常微分方程数值解
  • 几乎可以解决所有常微分方程组数值解问题。 论坛中不少朋友执着于利用MATLAB自置函数solve求解微分方程。其实solve函数存在局限性 (比如系统输入是时变的,solve函数就不易处理),也不利于大家掌握算法。 ...
  • 最近同学毕设需要求解循坏摆的微分方程,我在帮忙过程中学习了一下常微分方程的解析解和数值解的求法,在此分享。以下讲解遵循Matlab官方文档提供的方程和写法。(强烈建议大家有问题多看官方文档,非常有用)介绍一下...
  • 微分方程数值解系列博文: 偏微分方程数值解(一):定解问题 & 差分解法 偏微分方程数值解(二): 一维状态空间的偏微分方程的 MATLAB 解法 偏微分方程数值解(三): 化工应用实例 ----------触煤反应...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 3,240
精华内容 1,296
关键字:

常微分方程组的数值解