zoukankan      html  css  js  c++  java
  • MATLAB中ode23函数,龙格库塔函数

     今天说一说MATLAB中ode23函数的原理,在网上看了好多,但是不知道是怎么计算的,就知道是那么用的,但是最后结果咋回事不知道,今天来讲一讲是怎么计算的。

    首先来个程序:

    function f=eg6fun(t,y) 
    f=-y^3-2;  
    end

    上面是我定义的一个函数,看着挺简单的哈!不多说了。

    [t,y]=ode23(@eg6fun,[0,30],1);

    这句话是我用ode23调用的语句,先说一下,这里eg6fun是我上面函数的名称,也就是ode23主要计算的就是这个函数的微分。[0,30]为t的范围,这里t没有什么太大的作用,只是为了计算步长用的,之后我把运行后的t和y数据粘在这里,我们发现,在MATLAB中步长并不是固定的。这里应该是用一个什么函数求得,我没查,感兴趣的自己查一下。1为y的初值,也就是我们常常说的y0。

    先粘上实验结果,我们在分析怎么来的:

    下面的是t的值,这里MATLAB将t在[0,30]区间分成了67份,我这里只粘了一部分:

    0
    0.0266666666666667
    0.0974376132058027
    0.178185989598692
    0.270160382755732
    下面是y的结果,y最后也是一个[1,67]的矩阵:

    1
    0.922959859735161
    0.740501273051361
    0.556672216994644
    0.363414446549133

    下面我们来说是怎么计算的吧!看下面的图,这个是我在数值分析书上照的,其实ode23就是龙格库塔函数的应用,而龙格库塔函数就是根据欧拉法得来的,看下图:

    上面图片中有三个公式,第一个公式h后面括号中的内容就是要求积分的函数,就是我们程序中的eg6fun。那么就好办了,把图中公式中的括号中的内容换成我们的公式也就是

    -y^3-2

    然后计算就好了。这里h为步长,也就是我们程序中t的步长,我们可以看到第一次t为0.0266666666666667,而下一次的步长为0.0974376132058027-0.0266666666666667,只要这么一步一步计算就好了。

    (这里看图中黑色笔手写的公式)

    这里计算一步来表示计算的大概过程:

    例如: (1)计算Yp=y1+h * (-y1^3 - 2) = 1 - 0.027*3 = 0.919

       (2)       Yc=y1+h * (-Yp^3 - 2) = 1 - 0.027*(0.919^3 - 2) = 0.925

       (3)        Y(n+1) = 1/2 * (Yp + Yc) = 1/2 * (0.919 + 0.925) = 0.922

      因为这里我们保留精度为3位小数,可能计算的有些误差。还有一点需要注意,龙格库塔函数是对欧拉方法进行的改进,其实龙格库塔函数的精度要比欧拉方法更高。因此这里计算有些许误差。但是大概的过程就是这样的。

              上面的内容是之前写的,讲解的是欧拉算法计算微分的过程,其实龙格-库塔方法后来在书中看到,下面介绍一下龙格库塔方法:

           MATLAB中的ode23就是用的二阶的龙格库塔方法,就是图中3.6的三个公式,这里h为步长,上面给出的t,c1和c2是系数,这个系数取值不是固定的,MATLAB中是啥我也不是确定,但是书中最后给的是c1=0,c2=1,λ2和μ21取值1/2。这样一来,计算一波:y1=1;求y2,将y1带入公式中的yn,这里没有x,所以有x的项可以忽略

    k1=-3;

    k2=f(1-(1/2)*0.0267*3)=f(0.96)=-2.88

    y2=1-0.0267*2.88=0.923

           y2求出,其余的过程都是这样求得。ode45是四阶龙格库塔函数,下图为4阶求法,这里不再做介绍:

           到此MATLAB中ode23的计算方法已经讲解完了,当然,ode45跟这个应该类似,就是ode45比ode23更精确一点,在MATLAB中,如果我们用ode45会发现,t在[0,30]间分成了167份,很明显精度提高了。其实MATLAB中有好多的函数都是用到了数值分析中的内容,而数值分析就是用我们的笨方法来计算数值的一种工具,这是我自己定义的哈,通过减少误差来使计算出来的数据更准确。

  • 相关阅读:
    [Erlang05]gen_server怎么去写eunit?
    [Erlang04]为什么有了rpc还有net_kernel:connect/1?
    [Erlang03]Erlang有哪些好用的静态分析工具?
    [SIP01]SIP Header Fields里面各字段用途
    [SIP00]SIP 概念总结
    [Erlang00]:gen_server:reply/2
    Makefile教程
    Linux 环境下开发 STM32
    Ubuntu 18.04 + ROS Melodic + TurtleBot3仿真
    Ubuntu系统鼠标不能点击
  • 原文地址:https://www.cnblogs.com/Tang-YH/p/14327520.html
Copyright © 2011-2022 走看看