精华内容
下载资源
问答
  • 1.dx/dt=cos(y),t∈[0,T] dy/dt=sin(x),t∈[0,T] x(0)=0,y(0)=0 2.用梯形公式,显式欧拉,隐式欧拉,改进欧拉求解 3.有注释再好不过 4.可以有偿
  • 为了研究一阶线性模糊微分方程组的模糊初值问题,提出了模糊微分方程组的刻画方程组和关联解的概念,讨论了精确初值对刻画方程组解的影响,利用精确初值与关联解之间的关系,定义了模糊微分方程组初值问题解,同时给出了...
  • 常微分方程组数值求解;主要内容 数值求解常微分方程组函数概述 非刚性/刚性常微分方程问题求解 隐式微分方程组求解 微分代数方程(DAE)与延迟微分方程(DDE)求解 边值问题求解 ;第一节数值求解常微分方程组函数概述;一...
  • 简单的常微分方程()初值问题

    千次阅读 2015-08-17 19:03:42
    问题A,B,C,D,EA,B,C,D,E 五个化合物,AA 反应后生成 BB,BB 反应后生成CC,CC 反应后生成 DD,DD 反应后生成 EE,即A—→B—→C—→D—→EA—\rightarrow B—\rightarrow C—\rightarrow D—\rightarrow E ...

    问题

    ABCDE 五个化合物,A 反应后生成 BB 反应后生成CC 反应后生成 DD 反应后生成 E,即

    ABCDE

    反应动力学表明都是一级反应,已知 k1k2k3k4 是相应的反应速率常数。ABCDE 的浓度随时间变化的反应速率方程如下:
    1.2.3.4.5.dCAdt=k1CAdCBdt=k1CAk2CBdCCdt=k2CBk3CCdCDdt=k3CCk4CDdCEdt=k4CD

    当时间 t=0 时,A 的初始浓度为 A0,而 BCDE 的初始浓度都为 0

    A,B,C,D,E 浓度随时间变化的函数。

    解答

    简单的常微分方程(组)初值问题。可以顺次求解,也可以统一成方程组求。因为计算比较繁琐,用Maple或mathematica计算比较方便。在微分方程符号解方面,总体上看Maple比Mathematica胜出一筹。但在当前这个问题上,两种软件求解都易如反掌。

    Mathematica顺次求解ODE初值问题的代码如下:

    ClearAll["Global`*"]
    a[t_]:=Evaluate[a[t]/.(DSolve[{a'[t]==-k1 a[t],a[0]==A0},a,t]//FullSimplify)[[1,1]]]
    b[t_]:=Evaluate[b[t]/.(DSolve[{b'[t]==k1*a[t]-k2 b[t],b[0]==0},b,t]//FullSimplify)[[1,1]]]
    c[t_]:=Evaluate[c[t]/.(DSolve[{c'[t]==k2*b[t]-k3 c[t],c[0]==0},c,t]//FullSimplify)[[1,1]]]
    d[t_]:=Evaluate[d[t]/.(DSolve[{d'[t]==k3*c[t]-k4 d[t],d[0]==0},d,t]//FullSimplify)[[1,1]]]
    e[t_]:=Evaluate[e[t]/.(DSolve[{e'[t]==k4*d[t],e[0]==0},e,t]//FullSimplify)[[1,1]]]

    得到结果:

    CA(t)=CB(t)=CC(t)=CD(t)=CE(t)=A0ek1tA0k1ek1t(e(k1k2)t1)k1k2A0k1k2(ek3t(k1k2)ek2t(k1k3)+ek1t(k2k3))(k1k2)(k1k3)(k2k3)A0k1k2k3(ek4t(k1k4)(k4k2)(k4k3)+ek3t(k1k3)(k2k3)(k4k3)+ek2t(k1k2)(k3k2)(k4k2)+ek1t(k1k2)(k1k3)(k4k1))A0(k2k3k4ek1t(k1k2)(k1k3)(k1k4)k1k3k4ek2t(k1k2)(k2k3)(k2k4)k1k2k4ek3t(k1k3)(k3k2)(k3k4)k1k2k3ek4t(k1k4)(k4k2)(k4k3))

    最长的 CD(t) 只好贴出来:

    这里写图片描述

    CE(t):

    这里写图片描述

    这五个浓度之和是常数 A0 ,它们导数之和为 0 可以帮助验算。

    思考

    这里得到的解析解,乍一看似乎很方便。实际上,如果 ki=kj,(ij{1,2,3,4}) ,似乎出现了分母为 0 的情形。比如 k2=k3 时…… 难道是该动力学方程要求动力学常数不能相等吗?

    如何解决?从极限的角度,或者,直接通过微分方程特定条件下求解的方式

    展开全文
  • 我们采用了Duan-Rach-Wazwaz改进的Adomian分解方法来数值求解非线性常微分方程组的初值问题。 为了确认该方法的实用性,鲁棒性和可靠性,在某些情况下,我们将改进的Adomian分解方法与MATHEMATICA解决方案以及四阶...
  • 用 GSL 解常微分方程初值问题

    千次阅读 2012-05-12 12:22:46
    用 GSL 解常微分方程初值问题 GNU Scientific Library (GSL) 是一个用于科学计算...n 维的常微分方程组一般可以描述为: \[ \frac{{dy_i (t)}}{{dx}} = f_i (t,y_1 (t), \ldots y_n (t)) \] 其中 i = 1,\ldots,n

    用 GSL 解常微分方程初值问题

    GNU Scientific Library (GSL) 是一个用于科学计算的 C 语言类库。这里介绍如何用 GSL 来进行常微分方程初值问题的计算。

    n 维的常微分方程组一般可以描述为:
    \[
    \frac{{dy_i (t)}}{{dx}} = f_i (t,y_1 (t), \ldots y_n (t))
    \]
    其中 i = 1,\ldots,n
    相应的初值条件为 y_i (t_0) = Y_{i0}

    GSL 提供了许多的方法用来解决常微分方程的初值问题,包括底层的方法如 Runge-Kutta 法和 Bulirsch-Stoer 法,也包括高层的算法如自适应步长控制等。

    GSL 中用一个结构体 gsl_odeiv_system 来描述常微分方程组,

    typedef struct  
    {
      int (* function) (double t, const double y[], double dydt[], void * params);
      int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params);
      size_t dimension;
      void * params;
    }
    gsl_odeiv_system;

    function 是一个函数指针,指向用户定义的函数,t, y[], params作为输入参数,dydt 返回 n 维的向量场 f_i(t, y, params),计算成功时,函数应该返回 GSL_SUCCESS;

    jacobian 也是一个函数指针,指向用户定义的函数,t, y[], params作为输入参数,dfdt 返回 \partial f_i /\partial t,dfdy 返回 Jacobian 矩阵 J_{ij}, 存储方式为J(i,j) = dfdy[i * dimension + j],计算成功时,函数应该返回 GSL_SUCCESS;

    对于一些简单的算法并不需要 Jacobian 矩阵,此时可以将 jacobian 指向空(NULL).但是大多数高效的算法都要用到 Jacobian 矩阵,因此在能提供 Jacobian 矩阵时应该尽量的提供。

    dimension 描述系统的维数。

    params 为一个任意的指针,具体的作用在下面的例子中讲解。

    下面分为三个部分介绍 GSL 中提供的函数。

    1 步进函数

    步进函数是用户可以调用的最底层的功能函数,其他的高层函数也都是通过调用步进函数来完成最后的计算工作的。

    用来分配和释放所需空间

    gsl_odeiv_step * gsl_odeiv_step_alloc(const gsl_odeiv_step_type * T, size_t dim);
    int  gsl_odeiv_step_reset(gsl_odeiv_step * s);
    void gsl_odeiv_step_free(gsl_odeiv_step * s);

    其中 gsl_odeiv_step_type 可以为以下值:

    gsl_odeiv_step_rk2;    //二阶 Runge-Kutta 法
    gsl_odeiv_step_rk4;    //四阶经典 Runge-Kutta 法
    gsl_odeiv_step_rkf45;  //Runge-Kutta-Fehlberg 法
    gsl_odeiv_step_rkck;  //Runge-Kutta Cash-Karp
    gsl_odeiv_step_rk8pd; //Runge-Kutta Prince-Dormand
    gsl_odeiv_step_rk2imp; //隐式的 Runge-Kutta 
    gsl_odeiv_step_rk2simp;
    gsl_odeiv_step_rk4imp; //Implicit 4th order Runge-Kutta at Gaussian points
    gsl_odeiv_step_bsimp; //隐式的 Bulirsch-Stoer method of Bader and Deuflha (需要 Jacobian 矩阵)
    gsl_odeiv_step_gear1; //M=1 implicit Gear method
    gsl_odeiv_step_gear2; //M=2 implicit Gear method

    返回一些所用算法的信息

    const char * gsl_odeiv_step_name(const gsl_odeiv_step *);
    unsigned int gsl_odeiv_step_order(const gsl_odeiv_step * s);

    具体的步进计算函数如下

    int  gsl_odeiv_step_apply(gsl_odeiv_step *s, double t, double h, double y[], double yerr[], const double dydt_in[], double dydt_out[], const gsl_odeiv_system * dydt);

    计算 t -> t +h, t + h 时刻的值存储在 y[] 中,yerr[] 是对 y[] 计算误差的估计值。如果 dydt_in[] 不为 null ,dydt_in[] 应为 t 时刻的对时间的导数,如果没有提供 dydt_in[] 的话,这个函数的内部会作相应的计算。dydt_out[] 输出t + h 时刻对时间的导数。

    下面是一个简单的例子:
    考虑二阶的非线性Van der Pol 方程,
    \[
    x''(t) + \mu x'(t)(x(t)^2  - 1) + x(t) = 0
    \]
    这个方程可以被化简为如下的常微分方程组:

    y_1' = -y_0 + \mu y_1 (1 - {y_0}^2)
    y_0' = y_1

    Jacobian 矩阵为:
    \[
    J = \left( {\begin{array}{*{20}c}
       0 & 1  \\
       { - 1 - 2\mu y_0 y_1} & { \mu (1-{y_0}^2) }  \\
    \end{array}} \right)
    \]

    #include <stdio.h>
    #include <gsl/gsl_errno.h>
    #include <gsl/gsl_matrix.h>
    #include <gsl/gsl_odeiv.h>
    
    int func (double t, const double y[], double f[], void *params)
    {// 定义微分方程组
     double mu = *(double *)params;
     f[0] = y[1];
     f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1);
     return GSL_SUCCESS;
    }
    
    int jac (double t, const double y[], double *dfdy, double dfdt[], void *params)
    {
     double mu = *(double *)params;
     gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 2, 2);
     gsl_matrix * m = &dfdy_mat.matrix;
     gsl_matrix_set (m, 0, 0, 0.0);
     gsl_matrix_set (m, 0, 1, 1.0);
     gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0);
     gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0));
     dfdt[0] = 0.0;
     dfdt[1] = 0.0;
     return GSL_SUCCESS;
    }
    
    int main (void)
    {
     const gsl_odeiv_step_type * T = gsl_odeiv_step_rk8pd;
     gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 2);
    
     double mu = 10;
     gsl_odeiv_system sys = {func, jac, 2, &mu};
     double t = 0.0, t1 = 100.0;
     double h = 1e-6; // 步进为 h
     double y[2] = { 1.0, 0.0 }; //初值
     double y_err[2];
     double dydt_in[2], dydt_out[2];
     
     GSL_ODEIV_FN_EVAL(&sys, t, y, dydt_in); //计算初始时刻的dy/dt
     while (t < t1)
     {
      int status = gsl_odeiv_step_apply (
      s, t, h, y, y_err, dydt_in, dydt_out, &sys);
      if (status != GSL_SUCCESS) break;
      dydt_in[0] = dydt_out[0];
      dydt_in[1] = dydt_out[1];
      t += h;
      printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
     }
     gsl_odeiv_step_free (s);
     return 0;
    }

    程序很简单,不需要再解释了。
    这种方法的缺点是积分步长需要自己确定,步长太大计算结果会有很大的误差,步长太小效率低下,而且步长太小误差可能反而增大。

    因此,一般我们都会选用更高级的算法,也就是自适应步长算法。

    2 自适应步长控制

    gsl_odeiv_control * gsl_odeiv_control_standard_new(double eps_abs, double eps_rel, double a_y, double a_dydt);
    gsl_odeiv_control * gsl_odeiv_control_y_new(double eps_abs, double eps_rel);
    gsl_odeiv_control * gsl_odeiv_control_yp_new(double eps_abs, double eps_rel);
    gsl_odeiv_control * gsl_odeiv_control_scaled_new(double eps_abs, double eps_rel, double a_y, double a_dydt, const double scale_abs[], size_t dim);

    这几个函数大同小异,都是检查一个步进函数步进后的误差,与用户设定的误差水平进行比较,然后决定是否要调整步长。具体各个参数的含义请看手册。

    gsl_odeiv_control * gsl_odeiv_control_alloc(const gsl_odeiv_control_type * T);

    分配空间给新的步长控制函数,一般用户很少使用,除非自定义了新的步长控制函数

    int gsl_odeiv_control_init(gsl_odeiv_control * c, double eps_abs, double eps_rel, double a_y, double a_dydt);

    初始化步长控制函数

    void gsl_odeiv_control_free(gsl_odeiv_control * c);

    释放相应的空间

    int gsl_odeiv_control_hadjust (gsl_odeiv_control * c, gsl_odeiv_step * s, const double y0[], const double yerr[], const double dydt[], double * h);

    这个函数用来计算一个合适的步长 h, 步长增大时返回GSL_ODEIV_HADJ_INC,步长件小时返回GSL_ODEIV_HADJ_DEC,步长不做调整时返回GSL_ODEIV_HADJ_NIL。

    const char * gsl_odeiv_control_name(const gsl_odeiv_control * c);

    返回步长控制函数的名称。

    下面仍然以上面提到的 Van der Pol 方程为例。

    int main (void)
    {
     const gsl_odeiv_step_type * T = gsl_odeiv_step_rk8pd;
     gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 2);
     gsl_odeiv_control * c = gsl_odeiv_control_y_new (1e-5, 0.0);
     double mu = 10;
     gsl_odeiv_system sys = {func, jac, 2, &mu};
     double t = 0.0, t1 = 100.0;
     double h = 1e-2; // 步进为 h
     double *ph = &h;
     double y[2] = { 1.0, 0.0 }; //初值
     double y_const[2];
     double y_err[2];
     double dydt_in[2], dydt_out[2];
     int status;
     
     GSL_ODEIV_FN_EVAL(&sys, t, y, dydt_in); //计算初始时刻的dy/dt
     while (t < t1)
     {
      gsl_odeiv_step_apply (s, t, h, y_const, y_err, dydt_in, dydt_out, &sys);
      status = gsl_odeiv_control_hadjust (c, s, y_const, y_err, dydt_in, ph);
      gsl_odeiv_step_reset (s);
      
      status = gsl_odeiv_step_apply (s, t, h, y, y_err, dydt_in, dydt_out, &sys);
      
      y_const[0] = y[0];
      y_const[1] = y[1];
        
      if (status != GSL_SUCCESS) break;
      dydt_in[0] = dydt_out[0];
      dydt_in[1] = dydt_out[1];
      t += h;
      printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
     }
     gsl_odeiv_control_free(c);
     gsl_odeiv_step_free (s);
     return 0;
    }

    3 进化算法

    进化算法是最高层的算法,它结合了步进函数和步长控制函数。它可以自动调整步长。

    用来分配和释放所需空间

    gsl_odeiv_evolve * gsl_odeiv_evolve_alloc (size t dim);
    int gsl_odeiv_evolve_reset (gsl odeiv evolve * e);
    void gsl_odeiv_evolve_free (gsl odeiv evolve * e);
    int gsl_odeiv_evolve_apply (gsl odeiv evolve * e, gsl odeiv control * con, gsl odeiv step * step, const gsl odeiv system * dydt, double * t, double t1, double * h, double y[]);

    这个函数自动计算下一步的值,计算过程中可能会调整 h ,但是会保证 t + h <= t1,因此多次计算后会精确的到达t1。

    int main (void)
    {
     const gsl_odeiv_step_type * T = gsl_odeiv_step_rk8pd;
     gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 2);
     gsl_odeiv_control * c = gsl_odeiv_control_y_new (1e-6, 0.0);
     gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (2);
     double mu = 10;
     gsl_odeiv_system sys = {func, jac, 2, &mu};
     double t = 0.0, t1 = 100.0;
     double h = 1e-6;
     double y[2] = { 1.0, 0.0 };
     while (t < t1)
     {
      int status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t, t1, &h, y);
      if (status != GSL_SUCCESS) break;
      printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
     }
     gsl_odeiv_evolve_free (e);
     gsl_odeiv_control_free (c);
     gsl_odeiv_step_free (s);
     return 0;
    }

    经常我们需要微分方程的解在某些特定时刻的值(比如等时间间隔的一系列时刻的值),这时如果只是简单的采用进化算法是不行的,因为我们无法控制步长。下面的代码给出了一种解决办法。

    for (i = 1; i <= 100; i++)
    {
     double ti = i * t1 / 100.0;
     while (t < ti)
     {
      gsl_odeiv_evolve_apply (e, c, s, &sys, &t, ti, &h, y);
     }
     printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
    }




     

    展开全文
  • 常微分方程初值问题数值解法 一、问题提出 科学计算中经常遇到微分方程(初值问题,需要利用Euler法,改进Euler法,Rung-Kutta方法求其数值解,诸如以下问题: (1) (3) (4)利用四阶标准R- K...

     常微分方程初值问题数值解法

    一、问题提出

    科学计算中经常遇到微分方程(组)初值问题,需要利用Euler法,改进Euler法,Rung-Kutta方法求其数值解,诸如以下问题:

    (1)

    (3)

     

    (4)利用四阶标准R- K方法求二阶方程初值问题的数值解

    二、问题分析

        使用Euler法求解,运算程序简单,但是计算结果准确度不高。使用改进的Euler法求解过程相对复杂,但是准确度会更高。准确度最高的是四阶龙格库塔法,求解步骤也是最复杂的。问题(1)使用Euler求解,并与准确解对比。问题(3)使用改进的Euler法求解。问题(4)(I)(IV)使用四届标准龙格库塔法求解。

    三、实验程序及注释

    1.Euler法:

    View Code
    //differential.cpp
    #include <iostream>
    using namespace std;
    #include <math.h>

    float h, //步长
    x00, //x的取值范围
    xn, // x0< x < xn
    y00; //初值
    float x[21];//结果
    float y[21];

    void Euler()
    {
    init();
    x[0] = x00;
    y[0] = y00;
    for (int i=1; i<=(xn-x00)/h+1; i++)
    {
    x[i] = x[i-1] + h;
    y[i] = y[i-1] + h * f(x[i-1],y[i-1]);
    }
    }



    2.解问题(1):

    View Code
    //常微分方程y’= f(xn,yn)
    float f(float x, float y)
    {
    return 4*x/y-x*y;
    //return y-2*x/y;//课本上的例题,由于检测程序的正确性
    }
    //初始化初值和x的范围
    void init()
    {
    x00=0;
    xn=2;
    y00=3;
    }
    //求常微分方程的标准解
    float pf(float x)
    {
    return sqrt(4+5*exp(-x*x));
    }
    int main()
    {
    h=0.1;
    Euler();
    cout << "步长 h = "
    << h << endl;
    for (int i=1; i<=(xn-x00)/h+1; i++)
    {
    cout << "x" << i << "=" << x[i] << "\t"
    << "y" << i << "=" << y[i] << "\t"
    << "y("<< x[i] << ")=" << pf(x[i]) << endl;
    }
    h=0.2;
    Euler();
    cout << "步长 h = "
    << h << endl;
    for (i=1; i<=(xn-x00)/h+1; i++)
    {
    cout << "x" << i << "=" << x[i] << "\t"
    << "y" << i << "=" << y[i] << "\t"
    << "y("<< x[i] << ")=" << pf(x[i]) << endl;
    }
    h=0.4;
    Euler();
    cout << "步长 h = "
    << h << endl;
    for (i=1; i<=(xn-x00)/h+1; i++)
    {
    cout << "x" << i << "=" << x[i] << "\t"
    << "y" << i << "=" << y[i] << "\t"
    << "y("<< x[i] << ")=" << pf(x[i]) << endl;
    }
    return 0;
    }



    3.改进的Euler法:

    View Code
    #include <iostream>
    #include <iomanip>
    using namespace std;

    float h, //步长
    x00, //x的取值范围
    xn, // x0< x < xn
    y00; //初值
    float y10,y20,y30;//初值

    float x[21];//结果
    float y1[21],y2[21],y3[21];
    //改进的Euler法
    void newEuler()
    {
    init();
    x[0] = x00;
    y1[0] = y10;
    y2[0] = y20;
    y3[0] = y30;

    for (int i=1; i<=(xn-x00)/h+1; i++)
    {
    x[i] = x[i-1] + h;
    float yp=y1[i-1]+h*f1(x[i-1],y1[i-1],y2[i-1],y3[i-1]);
    float yc=y1[i-1]+h*f1(x[i],yp,y2[i-1],y3[i-1]);
    y1[i] = 0.5*(yp+yc);
    yp=y2[i-1]+h*f2(x[i-1],y1[i-1],y2[i-1],y3[i-1]);
    yc=y2[i-1]+h*f2(x[i],y1[i-1],yp,y3[i-1]);
    y2[i] = 0.5*(yp+yc);
    yp=y3[i-1]+h*f3(x[i-1],y1[i-1],y2[i-1],y3[i-1]);
    yc=y3[i-1]+h*f3(x[i],y1[i-1],y2[i-1],yp);
    y3[i] = 0.5*(yp+yc);
    }
    }



    4.解问题(3):

    View Code
    float f1(float x, float y1, float y2, float y3)
    {
    return y2;
    }
    float f2(float x, float y1, float y2, float y3)
    {
    return -y1;
    }
    float f3(float x, float y1, float y2, float y3)
    {
    return -y3;
    }
    void init()
    {
    x00=0;
    xn=0.15;
    y10=-1;
    y20=0;
    y30=1;
    }
    int main()
    {
    h=0.01;
    newEuler();
    cout << "步长 h = "
    << h << endl;
    cout.precision(6);
    for (int i=0; i<=(xn-x00)/h; i+=5)
    {
    cout << "x(" << setw(2) << i << ")=" << x[i] << ""
    << "y1(" << setw(2) << x[i] << ")=" << y1[i] << ""
    << "y2(" << setw(2) << x[i] << ")=" << y2[i] << ""
    << "y3(" << setw(2) << x[i] << ")=" << y3[i] << endl;
    }
    return 0;
    }



    5.龙格库塔法:

    View Code
    #include <iostream>
    #include <iomanip>
    using namespace std;

    float h, //步长
    x00, //x的取值范围
    xn; // x0< x < xn
    float y10;
    float y20;//初值
    float x[51];//结果
    float y1[51],y2[51];
    //龙格库塔法
    void L_k()
    {
    float L[5];
    init();
    x[0]=x00;
    y1[0]=y10;
    y2[0]=y20;
    for (int i=1; i<=(xn-x00)/h; i++)
    {
    x[i] = x[i-1] + h;
    L[1] = f2(x[i-1],y1[i-1],y2[i-1]);
    L[2] = f2(x[i-1]+h/2,y1[i-1]+h*y2[i-1]/2,y2[i-1]+h*L[1]/2);
    L[3] = f2(x[i-1]+h/2,y1[i-1]+h*y2[i-1]/2+h*h*L[1]/4,y2[i-1]+h*L[2]/2);
    L[4] = f2(x[i-1]+h,y1[i-1]+h*y2[i-1]+h*h*L[2]/2,y2[i-1]+h*L[3]);
    y1[i] = y1[i-1] + h*y2[i-1] + h*h*(L[1]+L[2]+L[3])/6;
    y2[i] = y2[i-1] + h*(L[1]+2*L[2]+2*L[3]+L[4])/6;
    }
    }



    6.解问题(4)(I):

    View Code
    //  y1'=y2
    // y2'=3*y2-*y1
    float f1(float x, float y1, float y2)
    {
    return y2;
    }
    float f2(float x, float y1, float y2)
    {
    return 3*y2-2*y1;
    }

    void init()
    {
    x00=0;
    xn=1;
    y10=0;
    y20=1;
    }

    int main()
    {
    h=0.02;
    L_k();
    cout << "步长 h = "
    << h << endl;
    cout.precision(6);
    for (int i=0; i<=(xn-x00)/h; i++)
    {
    cout << "x(" << setw(2) << i << ")=" << x[i] << ""
    << "y1(" << setw(2) << x[i] << ")=" << y1[i] << ""
    << "y2(" << setw(2) << x[i] << ")=" << y2[i] << endl;
    }
    return 0;
    }



    7.解问题(4)(IV):

    View Code
    //  y1'=y2
    // y2'=-sin(y1)
    float f1(float x, float y1, float y2)
    {
    return y2;
    }
    float f2(float x, float y1, float y2)
    {
    return -sin(y1);
    }

    void init()
    {
    x00=0;
    xn=4;
    y10=1;
    y20=0;
    }
    int main()
    {
    h=0.2;
    L_k();
    ……
    ……
    return 0;
    }



    四、实验数据结果及分析

    1.解问题(1):

         

    Euler法

    步长 h = 0.1

    x1=0.1 y1=3 y(0.1)=2.9917

    x2=0.2 y2=2.98333 y(0.2)=2.96714

    x3=0.3 y3=2.95048 y(0.3)=2.9274

    x4=0.4 y4=2.90264 y(0.4)=2.87415

    x5=0.5 y5=2.84166 y(0.5)=2.80963

    x6=0.6 y6=2.76995 y(0.6)=2.73649

    x7=0.7 y7=2.6904 y(0.7)=2.65766

    x8=0.8 y8=2.60615 y(0.8)=2.57613

    x9=0.9 y9=2.52044 y(0.9)=2.49485

    x10=1 y10=2.43643 y(1)=2.41648

    x11=1.1 y11=2.35697 y(1.1)=2.34329

    x12=1.2 y12=2.28438 y(1.2)=2.27698

    x13=1.3 y13=2.22038 y(1.3)=2.21869

    x14=1.4 y14=2.16592 y(1.4)=2.16894

    x15=1.5 y15=2.12124 y(1.5)=2.12767

    x16=1.6 y16=2.08591 y(1.6)=2.0944

    x17=1.7 y17=2.05898 y(1.7)=2.0683

    x18=1.8 y18=2.03922 y(1.8)=2.04837

    x19=1.9 y19=2.02523 y(1.9)=2.03353

    x20=2 y20=2.01571 y(2)=2.02276

    步长 h = 0.2

    x1=0.2 y1=3 y(0.2)=2.96714

    x2=0.4 y2=2.93333 y(0.4)=2.87415

    x3=0.6 y3=2.80776 y(0.6)=2.73649

    x4=0.8 y4=2.64178 y(0.8)=2.57613

    x5=1 y5=2.46136 y(1)=2.41648

    x6=1.2 y6=2.29411 y(1.2)=2.27698

    x7=1.4 y7=2.16199 y(1.4)=2.16894

    x8=1.6 y8=2.07467 y(1.6)=2.0944

    x9=1.8 y9=2.02774 y(1.8)=2.04837

    x10=2 y10=2.0079 y(2)=2.02276

    步长 h = 0.4

    x1=0.4 y1=3 y(0.4)=2.87415

    x2=0.8 y2=2.73333 y(0.8)=2.57613

    x3=1.2 y3=2.32696 y(1.2)=2.27698

    x4=1.6 y4=2.03513 y(1.6)=2.0944

    x5=2 y5=1.99055 y(2)=2.02276

     

        2.解问题(3):

    步长 h = 0.01

    x( 0)=0 y1( 0)=-1 y2( 0)=0 y3( 0)=1

    x( 5)=0.05 y1(0.05)=-0.999 y2(0.05)=0.04999 y3(0.05)=0.95123

    x(10)=0.1 y1(0.1)=-0.995502 y2(0.1)=0.09988 y3(0.1)=0.904839

    x(15)=0.15 y1(0.15)=-0.989514 y2(0.15)=0.149545 y3(0.15)=0.86071

    参考结果:

    对比看,结果相差不大,参考结果应该是使用四阶龙格库塔方法求解的。

    3.解问题(4)(I):

    步长 h = 0.02

    x( 0)=0 y1( 0)=0 y2( 0)=1

    x( 1)=0.02 y1(0.02)=0.0206094 y2(0.02)=1.06142

    x( 2)=0.04 y1(0.04)=0.0424763 y2(0.04)=1.12576

    x( 3)=0.06 y1(0.06)=0.0656603 y2(0.06)=1.19316

    x( 4)=0.08 y1(0.08)=0.0902238 y2(0.08)=1.26373

    x( 5)=0.1 y1(0.1)=0.116232 y2(0.1)=1.33763

    x( 6)=0.12 y1(0.12)=0.143752 y2(0.12)=1.415

    x( 7)=0.14 y1(0.14)=0.172856 y2(0.14)=1.49599

    x( 8)=0.16 y1(0.16)=0.203617 y2(0.16)=1.58074

    x( 9)=0.18 y1(0.18)=0.236112 y2(0.18)=1.66944

    x(10)=0.2 y1(0.2)=0.270422 y2(0.2)=1.76225

    x(11)=0.22 y1(0.22)=0.30663 y2(0.22)=1.85934

    x(12)=0.24 y1(0.24)=0.344825 y2(0.24)=1.9609

    x(13)=0.26 y1(0.26)=0.385097 y2(0.26)=2.06712

    x(14)=0.28 y1(0.28)=0.427543 y2(0.28)=2.17821

    x(15)=0.3 y1(0.3)=0.47226 y2(0.3)=2.29438

    x(16)=0.32 y1(0.32)=0.519353 y2(0.32)=2.41583

    x(17)=0.34 y1(0.34)=0.56893 y2(0.34)=2.54281

    x(18)=0.36 y1(0.36)=0.621104 y2(0.36)=2.67554

    x(19)=0.38 y1(0.38)=0.675991 y2(0.38)=2.81427

    x(20)=0.4 y1(0.4)=0.733716 y2(0.4)=2.95926

    x(21)=0.42 y1(0.42)=0.794405 y2(0.42)=3.11077

    x(22)=0.44 y1(0.44)=0.858192 y2(0.44)=3.26909

    x(23)=0.46 y1(0.46)=0.925216 y2(0.46)=3.43451

    x(24)=0.48 y1(0.48)=0.995622 y2(0.48)=3.60732

    x(25)=0.5 y1(0.5)=1.06956 y2(0.5)=3.78784

    x(26)=0.52 y1(0.52)=1.14719 y2(0.52)=3.97641

    x(27)=0.54 y1(0.54)=1.22867 y2(0.54)=4.17335

    x(28)=0.56 y1(0.56)=1.31418 y2(0.56)=4.37903

    x(29)=0.58 y1(0.58)=1.40389 y2(0.58)=4.59383

    x(30)=0.6 y1(0.6)=1.498 y2(0.6)=4.81811

    x(31)=0.62 y1(0.62)=1.59668 y2(0.62)=5.0523

    x(32)=0.64 y1(0.64)=1.70016 y2(0.64)=5.2968

    x(33)=0.66 y1(0.66)=1.80863 y2(0.66)=5.55205

    x(34)=0.68 y1(0.68)=1.92232 y2(0.68)=5.81851

    x(35)=0.7 y1(0.7)=2.04145 y2(0.7)=6.09665

    x(36)=0.72 y1(0.72)=2.16626 y2(0.72)=6.38696

    x(37)=0.74 y1(0.74)=2.29701 y2(0.74)=6.68996

    x(38)=0.76 y1(0.76)=2.43395 y2(0.76)=7.00617

    x(39)=0.78 y1(0.78)=2.57735 y2(0.78)=7.33617

    x(40)=0.8 y1(0.8)=2.72749 y2(0.8)=7.68052

    x(41)=0.82 y1(0.82)=2.88467 y2(0.82)=8.03984

    x(42)=0.84 y1(0.84)=3.04919 y2(0.84)=8.41474

    x(43)=0.86 y1(0.86)=3.22137 y2(0.86)=8.8059

    x(44)=0.88 y1(0.88)=3.40154 y2(0.88)=9.21397

    x(45)=0.9 y1(0.9)=3.59004 y2(0.9)=9.63969

    x(46)=0.92 y1(0.92)=3.78725 y2(0.92)=10.0838

    x(47)=0.94 y1(0.94)=3.99352 y2(0.94)=10.547

    x(48)=0.96 y1(0.96)=4.20926 y2(0.96)=11.0302

    x(49)=0.98 y1(0.98)=4.43487 y2(0.98)=11.5342

    x(50)=1 y1( 1)=4.67077 y2( 1)=12.0598



    3.解问题(4)(IV):

    步长 h = 0.2

    x( 0)=0 y1( 0)=1 y2( 0)=0

    x( 1)=0.2 y1(0.2)=0.983201 y2(0.2)=-0.167682

    x( 2)=0.4 y1(0.4)=0.933176 y2(0.4)=-0.331607

    x( 3)=0.6 y1(0.6)=0.851087 y2(0.6)=-0.48757

    x( 4)=0.8 y1(0.8)=0.739008 y2(0.8)=-0.630607

    x( 5)=1 y1( 1)=0.60009 y2( 1)=-0.754957

    x( 6)=1.2 y1(1.2)=0.438688 y2(1.2)=-0.854407

    x( 7)=1.4 y1(1.4)=0.260389 y2(1.4)=-0.923023

    x( 8)=1.6 y1(1.6)=0.0718545 y2(1.6)=-0.956155

    x( 9)=1.8 y1(1.8)=-0.119534 y2(1.8)=-0.951378

    x(10)=2 y1( 2)=-0.306182 y2( 2)=-0.90905

    x(11)=2.2 y1(2.2)=-0.480844 y2(2.2)=-0.832224

    x(12)=2.4 y1(2.4)=-0.637103 y2(2.4)=-0.725968

    x(13)=2.6 y1(2.6)=-0.769673 y2(2.6)=-0.596373

    x(14)=2.8 y1(2.8)=-0.874507 y2(2.8)=-0.449594

    x(15)=3 y1( 3)=-0.948738 y2( 3)=-0.291209

    x(16)=3.2 y1(3.2)=-0.990535 y2(3.2)=-0.125976

    x(17)=3.4 y1(3.4)=-0.998943 y2(3.4)=0.0420493

    x(18)=3.6 y1(3.6)=-0.973777 y2(3.6)=0.209154

    x(19)=3.8 y1(3.8)=-0.915597 y2(3.8)=0.371509

    x(20)=4 y1( 4)=-0.825779 y2( 4)=0.524739



    五、实验结论

       使用c语言写代码虽然比Matlab的代码长,但是还是c语言适合用来做练习,因为用c写更有结构。Matlab隐藏了许多细节部分,不适合对方法的深入理解,但是写代码速度快,代码更简练。

    转载于:https://www.cnblogs.com/hanxi/archive/2011/12/02/2272597.html

    展开全文
  • 我们重点是开发和实施新两步混合方法,用于直接求解一般二阶常微分方程。 在该方法开发中,采用幂级数作为基础函数,并且将出现的方程组微分系统并置在所有电网和离网点。 所得方程在选定点处插值。 然后,...
  • #先上代码后补笔记# #可以直接复制粘贴调用MATLAB函数代码!# 1. 朗格-库塔(Runge-Kutta)方法族 目前只实现了四阶Runge-Kutta方法。...% 朗格-库塔(Runge-Kutta)算法解常微分方程); % 输入四...

    #先上代码后补笔记#

    #可以直接复制粘贴调用的MATLAB函数代码!#

    1. 朗格-库塔(Runge-Kutta)方法族

    目前只实现了四阶Runge-Kutta方法。

    function [ YMat ] = Runge_Kutta( func, tvec, y_init, order )
    %   朗格-库塔(Runge-Kutta)算法解常微分方程(组);
    %   输入四个参数:函数句柄func(接收列向量、返回列向量),积分时间列向量tvec,初值行向量y_init,阶数order;
    %   输出一个参数:数值解,每一行对应积分时间列向量的一行,各列为变量一个分量。
    %   暂时只实现了阶数为4的(即order='4')
    col = size(y_init, 2); row = size(tvec, 1); YMat = zeros(row, col);
    YMat(1, :) = y_init;
    switch order
        case '4'
            for i=1:row-1
                y = YMat(i, :).'; h = tvec(i+1) - tvec(i); t = tvec(i);
                k1 = func(t, y);
                k2 = func(t + h/2, y + k1*h/2);
                k3 = func(t + h/2, y + k2*h/2);
                k4 = func(t + h, y + k3*h);
                YMat(i + 1, :) = (y + (k1 + 2*k2 + 2*k3 + k4)*h/6).';
            end
    end
    

      

    转载于:https://www.cnblogs.com/gentle-min-601/p/9638170.html

    展开全文
  • 2.3 常微分方程组的数值解法 知识要点 常微分方程初值问题---ode45,0de23 常微分方程边值问题---bvp4c 微分方程在化工模型中的应用 间歇反应器的计算 活塞流反应器的计算 全混流反应器的动态模拟 定态一维热传导问题...
  • 第9章-常微分方程初值问题数值解法9.1 引言9.2 简单数值方法9.3 龙格-库塔方法9.4 单步法收敛性与稳定性9.5 线性多步法9.6 线性多步法收敛性与稳定性9.7 一阶方程组与刚性方程组 9.1 引言 9.2 简单数值方法 ...
  • 1求解常微分方程MATLAB程序为 2求解常微分方程MATLAB程序为 3求常微分方程组通解MATLAB程序为 4求常微分方程组通解MATLAB程序为 5求解常微分方程MATLAB程序如下 6求解常微分方程解并画出解图形 ?...
  • 数值方法: 我们实验目标是解常微分方程其中包括几类问题一阶常微分初值问题高阶 常微分初值问题常微分方程组初值问题二阶常微分方程边值问题二阶线性常 微分方程边值问题 对待上面几类问题我们分别使用不同...
  • 输入所求微分方程组的初值 t0,x1,y1,x2,y2,x3,y3,vx1,vy1,vx2,vy2,vx3,vy3 输入所求微分方程组的微分区间[a,b] 输入所求微分方程组所分解子区间的个数step 输出文件可以直接用Matlab来作图 当时怎么写的我已经...
  • 如 可化为两个一阶微分方程 常微分方程组的解不完全由其方程决定,要确定方程的解还需要边界值条件。边界值条件分为初值问题和两点边值问题,本章主要讨论初值问题。 1. 龙格库塔法Runge-Kutta 通过...
  • 1 往期链接Chenglin Li:数值计算(四十五)Simulink求解非线性方程组​zhuanlan.zhihu.com2 算法介绍(如果公式格式不对,请刷新一下网址~)对于一般一阶常微分方程组,考虑初值问题 四阶Runge-Kutta方法可以表示...
  • 四阶龙格库塔方法求解一次常微分方程组一、写在前面二、四阶龙格库塔方法三、使用四阶龙格库塔方法求解一次常微分方程组 一、写在前面 龙格库塔方法是数值求解常微分非线性方程有利工具,计算精度较高,通过缩短...
  • 本节将介绍常微分方程初值问题的数值求解,主要内容分为三个部分:常微分方程数值求解概念、求解函数及刚性问题。 一、常微分方程数值求解一般概念 首先,凡含有参数,未知函数和未知函数导数 (或微分) 方程,...
  • 第五部分 常微分方程数值解;...常微分方程组odefile编写;高阶微分方程odefile编写;解算指令使用方法;常微分方程初值问题解算指令比较;ode解算指令选择(1;ode45和ode23比较;ode解算指令选择
  • 第五讲 常微分方程数值解 化工学院软件应用教科 2006-10 本章知识要点 数值计算 常微分方程初值问题 常微分方程边值问题 MATLAB 微分方程求解常微分方程的相关函数 ode45 ode23 bvp4c 微分方程在化工模型中应用 ...
  • 本文着重于开发一种具有块扩展的混合方法,用于直接求解一般三阶常微分方程的初值问题(IVP)。 幂级数用作IVP解决方案的基础函数。 将基函数的近似解插值到某些选定的离网点,同时将近似解的三阶导数并置在所有网格...
  • 一、问题的引入——兰彻斯特作战模型(作战双方的兵员数量)和军备竞赛模型(甲乙两...2. 初值问题方程组的特解 3. 齐次线性微分方程组解的结构(通解) 4. 非齐次线性微分方程组解的结构(通解+特解) ...
  • Matlab中ode系列函数的用法 Matlab的ode系列函数(统称odesolver)专门用于求解常微分方程的初值问题。ode是常微分方程(ordinary differential equation)的英文缩写。 Matlab中采用ode系列函数求解常微分...
  • 通过使用本文提出的方法(LOM)获得的结果表明,数值方法非常有效,适合求解分数阶常微分方程的初值和边值问题。 此外,还提供了一些数值示例,并对获得的结果与获得的分析结果进行了比较,证明了该方法的有效性。
  • PAGE / NUMPAGES 第7章 MATLAB解方程与函数极值 7.1 线性方程组求解 7.2 非线性方程数值求解 7.3 常微分方程初值问题的数值解法 7.4 函数极值 7.1 线性方程组求解 7.1.1 直接解法 1利用左除运算符直接解法 对于...
  • MATLAB 线性方程组求解 MATLAB 非线性方程数值求解 MATLAB 常微分方程初值问题的数值解法 MATLAB 最优化问题求解 6.1 线性方程组求解 6.1.1 直接解法 1 利用左除运算符直接解法 高斯消元法主元素消元法平方根法...
  • 1. 误差 2. 多项式插值与样条插值 3. 函数逼近 4. 数值积分与数值微分 5. 线性方程组的直接解法 6. 线性方程组的迭代解法 7. 非线性方程求根 8. 特征值和特征向量的计算 9. 常微分方程初值问题的数值解

空空如也

空空如也

1 2 3 4 5
收藏数 85
精华内容 34
关键字:

常微分方程组的初值问题