精华内容
下载资源
问答
  • 利用 MATLAB 求解常微分方程数值解 目录 1. 内容简介 把高等工程数学看了一遍增加对数学内容的了解对其中数值解法比较感兴趣 这大概是因为在其它各方面的学习和研究中经常会遇到数值解法的问题理解模型然后列出 微分...
  • 本文是自己写的关于怎样利用MATLAB求解常微分方程数值解的,文中从Euler法讲起,最后总结了常用的odeXX的用法及其原理,其中包含各个函数怎样使用的MATLAB代码
  • 常微分方程数值求解;主要内容 数值求解常微分方程组函数概述 非刚性/刚性常微分方程问题求解 隐式微分方程组求解 微分代数方程(DAE)与延迟微分方程(DDE)求解 边值问题求解 ;第一节数值求解常微分方程组函数概述;一...
  • 差分格式解决偏微分方程,用欧拉方法求解下列常微分方程的初值问题。 显示Euler格式,改进的Euler格式
  • 数值分析(数值计算) 数学建模 实验报告 matlab程序
  • 用Eular法解常微分方程组的数值解,使用了细胞数组,代码简洁,除注释外的有效代码只有二十行左右。(几年前上传的程序了,当时要20积分,现在为大家降到5个积分)
  • 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还支持带有时变项和额外参数的微分方程求解,这里不再赘述,大家可以自行参阅官方文档。

    展开全文
  • matlab常微分方程数值求解

    千次阅读 2020-09-01 21:08:13
    本节将介绍常微分方程初值问题的数值求解,主要内容分为三个部分:常微分方程数值求解的概念、求解函数及刚性问题。 一、常微分方程数值求解的一般概念 首先,凡含有参数,未知函数和未知函数导数 (或微分) 的方程,...

    本节将介绍常微分方程初值问题的数值求解,主要内容分为三个部分:常微分方程数值求解的概念、求解函数及刚性问题。
    一、常微分方程数值求解的一般概念
    首先,凡含有参数,未知函数和未知函数导数 (或微分) 的方程,称为微分方程,有时简称为方程,未知函数是一元函数的微分方程称作常微分方程,未知函数是多元函数的微分方程称作偏微分方程。微分方程中出现的未知函数最高阶导数的阶数,称为微分方程的阶。
    常微分方程数值求解,指研究求解初值问题各类数值方法的构造、理论分析与数值实现问题。研究的主要对象为一阶方程组初值问题:
    在这里插入图片描述
    其中,其中y=y(x)是未知函数,y(x0)=y0是初值条件,而f(x,y)是给定的二元函数。
    所谓其数值解法,就是求y(x)在离散结点xn处的函数近似值yn 的方法 ,yn≈y(xn)。这些近似值称为常微分方程初值问题的数值解。相邻两个结点之间的距离称为步长。
    这里主要介绍两种:单步法和多步法
    单步法:在计算y(n+1)时只用到前一步的y(n),因此在有了初值之后就可以逐步往下计算,其代表是龙格- - 库塔( Runge- - Kutta ) 法
    多步法:在计算y(n+1)时,除了用到前一步的值y(n)之外, , 还要用到y(n-p)( p=1,2 , … ,k,k>0)的值, , 即前面的k步。其代表就是亚当斯 (Adams) 法
    更多介绍可参见这个链接:
    https://wenku.baidu.com/view/cafe161b9a6648d7c1c708a1284ac850ad0204fc.html
    二、常微分方程求解函数
    MATLAB 提供了多个求常微分方程初值问题数值解的函数,一般调用格式为:
    [t,y]=solver(filename,tspan,y0,option)
    其中,t和y分别给出时间向量和相应的数值解。solver为求常微分方程数值解的函数。filename 是定义 f(t ,y) 的函数名,该函数必须返回一个列向量。
    tspan 形式为 [t0,tf],表示求解区间。 y0 是初始状态向量。Option 是可选参数,用于设置求解属性,常用的属性包括相对误差值 RelTol(默认值是10的-3次方)和绝对误差值 AbsTol( ( 默认值是10的-6次方)。
    常微分方程数值求解函数的统一命名格式:
    odennx
    ode是Ordinary Differential Equation 的缩写,是常微分方程的意思。
    nn 是数字,代表所用方法的阶数。例如,ode23采用2阶龙格- - 库塔( Runge- - Kutta )算法 ,用3阶公式做误差估计来调节步长,具有低等精度。ode45 采用4阶龙格- - 库塔算 法 ,用5阶公式做误差估计来调节步长,具有中等精度。
    x是字母,用于标注方法的专门特征。例如, ode15s 、ode23s 中的“s”代表( Stiff ),表示函数适用于刚性方程。
    下表列出了求常微分方程数值解的函数:

    求解函数采用方法适用场合
    ode232 阶或3阶 Runge- - Kutta 算法,低精度非刚性
    ode454 阶或5阶 Runge- - Kutta 算法,中精度非刚性
    ode113Adams 算法,精度可到10的-3次方至10的-6次方非刚性,计算时间比 ode45
    ode23t梯形算法适度刚性
    ode15sGear’s 反向数值微分算法,中精度刚性
    ode23s2阶 Rosebrock 算法,低精度刚性,当精度较低时,计算时间比 ode15s
    ode23tb梯形算法,低精度刚性,当精度较低时,计算时间比 ode15s

    先看一个简单例子,y‘=y+3x/x^2,初值y(0)=-2,求解区间为[1,4]。

    >> odefun=@(x,y) (y+3*x)/x^2;
    tspan=[1 4];
    y0=-2;
    [x y]=ode45(odefun,tspan,y0)
    plot(x,y)
    

    在这里插入图片描述 二、刚性问题
    有一类常微分方程,其解的分量有的变化很快,有的变化很慢,且相差悬殊,这就是所谓的刚性问题 (Stiff) 。
    对于刚性问题,数值解算法必须取很小步长才能获得满意的结果,导致计算量会大大增加。解决刚性问题需要有专门方法。非刚性算法可以求解刚性问题,只不过需要很长的计算时间。
    例、假如点燃一个火柴,火焰球迅速增大直至某个临界体积,然后,维持这一体积不变,原因是火焰球内部燃烧耗费的氧气和从球表面所获氧气达到平衡。其简化模型如下:
    y’=y2-y3,y(0)=L,0<=t<=2/L
    其中, y(t) 代表火焰球半径,初始半径是λ ,它很小。分析 λ 的大小对方程求解过程的影响。

    >>  L=0.01;
    f=@(t,y) y^2-y^3;
    tic;[ t,y ]=ode45(f,[0,2/L],L); toc
    disp (['ode45 计算的点数' num2str(length(t))]);
    时间已过 0.003704 秒。
    ode45 计算的点数157
    >>  L=1e-5;
    f=@(t,y) y^2-y^3;
    tic;[ t,y ]=ode45(f,[0,2/L],L); toc
    disp (['ode45 计算的点数' num2str(length(t))]);
    时间已过 1.010016 秒。
    ode45 计算的点数120565
    >> L=1e-5;
    f=@(t,y) y^2-y^3;
    tic;[ t,y ]=ode15s(f,[0,2/L],L); toc
    disp (['ode45 计算的点数' num2str(length(t))]);
    时间已过 0.189977 秒。
    ode45 计算的点数98
    

    注:tic 和 toc 函数 用来记录 微分方程求解 命令执行的时间 ,使用 tic 函数启动计时器,使用 toc 函数显示从计时器启动到当前所经历的时间。
    上述采用了三种不同方法,可以发现,第一种的运行结果表明这时常微分方程不算很刚性;第二种这时计算时间明显加长,计算的点数剧增,表明这时常微分方程表现为刚性;因此对于刚性问题,我们需要改变求解算法,我们选择以“s”结尾的函数,例如第三种方法他们专门用于求解刚性方程。计算时间大大缩短,计算的点数大大减少。
    常微分方程数值解法的研究已发展得相当成熟,理论上也颇为完善,小编由于能力有限,只能了解个大概,本文讲解也就写的是比较基础的一些方面,如果大家有更多需求,可以自己查阅相关资料。

    关于MATLAB的学习:

    大家可以关注我们的知乎专栏——数据可视化和数据分析中matlab的使用:
    https://zhuanlan.zhihu.com/c_1131568134137692160

    欢迎大家加入我们的MATLAB学习交流群:
    953314432
    扫码关注我们
    发现更多精彩
    在这里插入图片描述

    展开全文
  • MATLAB算法-求解微分方程数值解和解析解.ppt
  • 8.7 常微分方程数值解的若干Matlab函数文件; function [x,y]=euler(f,tspan,y0,n) % 解初值问题 =f(x,y,y(a)=y0 % 使用 n 步的 Euler 法步长 h=(b-a)/h a=tspan(1;b=tspan(2;h=(b-a)/h; x=(a+h:h:b; y(1)=y0+h*feval...
  • 主要介绍了python中sympy库求常微分方程的用法,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面随着小编来一起学习学习吧
  • 用ODE(ordinary differential equations)23和ODE45求解微分方程 ODE(ordinary differential equations)23 通过使用二阶和三阶隆哥-库塔方法积分微分方程得出

    用ODE(ordinary differential equations)23和ODE45求解微分方程

    ODE(ordinary differential equations)23

    [t,y]=ode23('func_name',[start_time,end_time],y(0))

    通过使用二阶和三阶隆哥-库塔方法积分微分方程得出

    function ydot=eq1(t,y)
    ydot=cos(t)

    [t,y]=ode23('eq1',[0 2*pi],2);----------就会产生一系列的数值解

    ydot =

         1


    ydot =

        0.9968


    ydot =

        0.9928


    ydot =

        0.9872


    ydot =

        0.9400


    ydot =

        0.9038


    ydot =

        0.8596


    ydot =

        0.7394


    ydot =

        0.6676


    ydot =

        0.5890


    ydot =

        0.3909


    ydot =

        0.2835


    ydot =

        0.1724


    ydot =

       -0.1176


    ydot =

       -0.2604


    ydot =

       -0.3977


    ydot =

       -0.6618


    ydot =

       -0.7709


    ydot =

       -0.8610


    ydot =

       -0.6018


    ydot =

       -0.6919


    ydot =

       -0.7723


    ydot =

       -0.8998


    ydot =

       -0.9450


    ydot =

       -0.9770


    ydot =

       -1.0000


    ydot =

       -0.9953


    ydot =

       -0.9799


    ydot =

       -0.9248


    ydot =

       -0.8846


    ydot =

       -0.8365


    ydot =

       -0.7234


    ydot =

       -0.6577


    ydot =

       -0.5865


    ydot =

       -0.4280


    ydot =

       -0.3429


    ydot =

       -0.2550


    ydot =

       -0.0510


    ydot =

        0.0524


    ydot =

        0.1553


    ydot =

        0.4529


    ydot =

        0.5868


    ydot =

        0.7063


    ydot =

        0.3516


    ydot =

        0.4448


    ydot =

        0.5334


    ydot =

        0.6932


    ydot =

        0.7628


    ydot =

        0.8245


    ydot =

        0.9143


    ydot =

        0.9477


    ydot =

        0.9730


    ydot =

        0.9932


    ydot =

        0.9983


    ydot =

         1

    >> [t,y]=ode23('eq1',[0 2*pi],2);

    ydot =

         1


    ydot =

        0.9968


    ydot =

        0.9928


    ydot =

        0.9872


    ydot =

        0.9400


    ydot =

        0.9038


    ydot =

        0.8596


    ydot =

        0.7394


    ydot =

        0.6676


    ydot =

        0.5890


    ydot =

        0.3909


    ydot =

        0.2835


    ydot =

        0.1724


    ydot =

       -0.1176


    ydot =

       -0.2604


    ydot =

       -0.3977


    ydot =

       -0.6618


    ydot =

       -0.7709


    ydot =

       -0.8610


    ydot =

       -0.6018


    ydot =

       -0.6919


    ydot =

       -0.7723


    ydot =

       -0.8998


    ydot =

       -0.9450


    ydot =

       -0.9770


    ydot =

       -1.0000


    ydot =

       -0.9953


    ydot =

       -0.9799


    ydot =

       -0.9248


    ydot =

       -0.8846


    ydot =

       -0.8365


    ydot =

       -0.7234


    ydot =

       -0.6577


    ydot =

       -0.5865


    ydot =

       -0.4280


    ydot =

       -0.3429


    ydot =

       -0.2550


    ydot =

       -0.0510


    ydot =

        0.0524


    ydot =

        0.1553


    ydot =

        0.4529


    ydot =

        0.5868


    ydot =

        0.7063


    ydot =

        0.3516


    ydot =

        0.4448


    ydot =

        0.5334


    ydot =

        0.6932


    ydot =

        0.7628


    ydot =

        0.8245


    ydot =

        0.9143


    ydot =

        0.9477


    ydot =

        0.9730


    ydot =

        0.9932


    ydot =

        0.9983


    ydot =

         1

    >> f=2+sin(t);
    >> plot(t,y,'o',t,f)

    这一系列数值解就是原函数的曲线上。y是数值解,f是解析解

    查看数值解与解析解的误差

     >> for i=1:1:size(y)
    err(i)=abs((f(i)-y(i))/f(i));
    end

    ODE(ordinary differential equations)45

    [t,w]=ode45('function',[start_time,end_time],y(0))

    使用高阶龙格-库塔公式

    >> [t,w]=ode45('eq1',[0 2*pi],2);

    ydot =

         1


    ydot =

        0.9968


    ydot =

        0.9927


    ydot =

        0.9488


    ydot =

        0.9369


    ydot =

        0.9203


    ydot =

        0.9203


    ydot =

        0.8640


    ydot =

        0.8307


    ydot =

        0.6180


    ydot =

        0.5732


    ydot =

        0.5146


    ydot =

        0.5146


    ydot =

        0.4031


    ydot =

        0.3449


    ydot =

        0.0379


    ydot =

       -0.0179


    ydot =

       -0.0876


    ydot =

       -0.0876


    ydot =

       -0.2118


    ydot =

       -0.2727


    ydot =

       -0.5567


    ydot =

       -0.6022


    ydot =

       -0.6564


    ydot =

       -0.6564


    ydot =

       -0.7458


    ydot =

       -0.7862


    ydot =

       -0.9387


    ydot =

       -0.9564


    ydot =

       -0.9745


    ydot =

       -0.9745


    ydot =

       -0.9949


    ydot =

       -0.9993


    ydot =

       -0.9621


    ydot =

       -0.9454


    ydot =

       -0.9203


    ydot =

       -0.9203


    ydot =

       -0.8640


    ydot =

       -0.8307


    ydot =

       -0.6180


    ydot =

       -0.5732


    ydot =

       -0.5146


    ydot =

       -0.5146


    ydot =

       -0.4031


    ydot =

       -0.3449


    ydot =

       -0.0379


    ydot =

        0.0179


    ydot =

        0.0876


    ydot =

        0.0876


    ydot =

        0.2118


    ydot =

        0.2727


    ydot =

        0.5567


    ydot =

        0.6022


    ydot =

        0.6564


    ydot =

        0.6564


    ydot =

        0.7458


    ydot =

        0.7862


    ydot =

        0.9387


    ydot =

        0.9564


    ydot =

        0.9745


    ydot =

        0.9745


    ydot =

        0.9836


    ydot =

        0.9875


    ydot =

        0.9990


    ydot =

        0.9997


    ydot =

         1


    ydot =

         1

    >> f=2+sin(t);
    >> plot(t,w,'0',t,f)

     


    ode45返回了45个数据点,而ode23返回7个数据点。


    求两个一阶方程组的数值解:

    \frac{\mathrm{dx} }{\mathrm{d} t}=-x^2+y

    \frac{\mathrm{dy} }{\mathrm{d} t}=-x+xy

    .m文件:

    function xdot=eqx(t,x)
    xdot=zeros(2,1)
    xdot(1)=-x(1)^2+x(2);
    xdot(2)=-x(1)-x(1)*x(2)

    >> [t,x]=ode45('eqx',[0 10],[0 1]);

    >> plot(t,x(:,1),t,x(:,2),'--')

     >> plot(x(:,1),x(:,2)),xlabel('x(1)'),ylabel('x(2)')

    求解二阶微分方程

    {y}''+16y=sin(4.3t)

    y(0)=y'(0)=0


    首先转换成一个一阶微分方程组:

    x_1=y,x_2=y'

    x_1'=y'=x_2;x_2'=y''=sin(4.3t)-16x_1


    .m文件

    function xdot=eqx(t,x)
    xdot=zeros(2,1)
    xdot(1)=x(2);
    xdot(2)=sin(4.3*t)-16*x(1);

    命令行窗口:

    >> [t,x]=ode45('eqx',[0 2*pi],[0,0]);

    >> plot(t,x(:,1),t,x(:,2),'--')

     >> plot(x(:,1),x(:,2)),xlabel('x(1)'),ylabel('x(2)')

     


    求解:

    x_1'=2x_1x_2;x_2'=-x_1^2;x_1(0)=x_2(0)=0


    .m文件:

    function xdot=eqx(t,x)
    xdot=zeros(2,1)----------------------------(分配存储空间)
    xdot(1)=2*x(1)*x(2);
    xdot(2)=x(1)^2;

    >>[t,x]=ode45('eqx',[0 2*pi],[0,0])

    >> plot(t,x(:,1),t,x(:,2),'--')

     

    >>[t,x]=ode45('eqx',[0 2*pi],[1,2])

    >> plot(t,x(:,1),t,x(:,2),'--')

     


    求解:

    y'=-2.3y;y(0)=0


    function ydot=eqx(t,y)
    ydot=-2.3*y;

    >> [t,x]=ode45('eqx',[0 2*pi],[1,2])

    >> plot(t,y)

     >> syms y(t);
    >> eqn=diff(y,t)+2.3*y==0;
    >> cond=(y(0)==0);
    >> s=dsolve(eqn,cond)
     
    s =
     
    0


    求解:

    y''-2y'+y=exp(-t);y(0)=2,y'(0)=0



    x_1=y;x_2=y'

    x_1'=y'=x_2

    x_2'=y''=exp(-t)-y+2y'=exp(-t)-x_1+2x_2


    .m文件:

    >> [t,x]=ode45('eqx',[0 2*pi],[2,0])

    >> plot(t,x(:,1),t,x(:,2),'--')

     >> syms y(t);
    >> eqn=diff(y,t,2)-2*diff(y,t)+y==exp(-t);
    >> Dy=diff(y,t);
    >> cond=[y(0)==2,Dy(0)==0];
    >> s=dsolve(eqn,cond)
     
    s =
     
    (exp(-t)*(7*exp(2*t) - 6*t*exp(2*t) + 1))/4

    >> ezplot(s)

     


    直接用dsolve:可以直接解出方程

    而用ode45可以解出一阶微分方程,或者二阶微分方程中状态变量的解,更加使用与线性系统理论中的状态空间表达式的各个状态变量随时间的变化情况。


     

    展开全文
  • *题目描述:*求常微分方程y’=y-2x/y && y(0)=1的数值解。解析解为:y=sqrt(1+2x) 编译环境vc++6.0: 代码如下: /* 2020151401 黄金 常微分方程数值解 1、y'=y-2x/y && y(0)=1; 2、解析解为y=...

    *题目描述:*求常微分方程y’=y-2x/y && y(0)=1的数值解。解析解为:y=sqrt(1+2x)
    编译环境vc++6.0
    代码如下:

    /*
    1、y'=y-2x/y && y(0)=1;
    2、解析解为y=sqrt(1+2x)
    */
    
    #include <stdio.h>
    #include <math.h>
    #define N 5                    //离散点个数,分别取5,8,10,对应步长0.2,0.125,0.1
    #define low 0					//定义域x区间下限
    #define high  1.0				//定义域x区间上限
    #define h  (high-low)/N        //步长
    
    /*************************4、解析解实现**********************/
    double AnalyticalSolution(double x) {
    	return sqrt(1 + 2 * x);
    }
    
    /*************************1、显示欧拉实现**********************/
    double EulerMethod(double x)
    {
    	//若x==0
    	if (fabs(x - low) < 1e-6)  //浮点数之间不能直接判断是否相等,则可以用 fabs(f1-f2)<=给定差值精度 来判断f1和f2是否相等
    		return 1;				//f(low)== f(0)==1
    	else return EulerMethod(x - h) + h * (EulerMethod(x - h) - (2 * (x - h) / EulerMethod(x - h)));
    }
    
    /*************************2、隐式欧拉实现**********************/
    //隐式欧拉实现由下面两个函数组成
    double ImplicitEulerMethod_1(double x)			/*隐示欧拉第一次迭代,
    用显示欧拉计算f(x)的初值,带入隐式欧拉方程 f(x)=f(x-1)+h*(f(x)-2*x/f(x))*/
    {
    	if (fabs(x - low) < 1e-6)	//浮点数之间不能直接判断是否相等,则可以用 fabs(f1-f2)<=给定差值精度 来判断f1和f2是否相等
    		return 1;
    	else
    		return ImplicitEulerMethod_1(x - h) + h * (EulerMethod(x) - 2 * x / EulerMethod(x));   //迭代
    }
    
    double ImplicitEulerMethod_2(double x)	/*隐示欧拉第2--n次迭代,
    用第一次迭代的结果f(x)作为初值带入隐式欧拉方程进行第二次迭代,
    第二次迭代的结果f(x)作为初值带入隐式欧拉方程进行第三次迭代.....
    直到本次迭代的结果y1与上次的迭代结果y0,满足|y1-y0|<1e-6,则认为收敛
    迭代结束
    */
    {
    	if (fabs(x - low) < 1e-6)  //浮点数之间不能直接判断是否相等,则可以用 fabs(f1-f2)<=给定差值精度 来判断f1和f2是否相等
    		return 1;
    	else {
    		double y0 = ImplicitEulerMethod_1(x), y1 = ImplicitEulerMethod_1(x);		//将第一次迭代的结果作为初值带入第二次迭代
    		do {
    			y0 = y1;
    			y1 = ImplicitEulerMethod_2(x - h) + h * (y0 - 2 * x / y0);
    		} while (fabs(y1 - y0) < 1e-6);		//本次迭代结果减去上一次迭代结果fabs(y1-y0)<=1e-6  认为是收敛
    		return y1;
    	}
    }
    /*************************3、两步欧拉实现**********************/
    double DoubleStepEuler(double x)			//两步欧拉
    {
    	if (fabs(x - low) < 1e-6)	//若x==0   第一个初值给定
    		return 1;
    	else if (fabs(x - low - h) < 1e-6)   //若x==low+h	 第二个初值用显示欧拉求出
    		return EulerMethod(x);
    	else return DoubleStepEuler(x - 2 * h) + 2 * h * (DoubleStepEuler(x - h) - 2 * (x - h) / DoubleStepEuler(x - h));
    }
    
    
    
    /*************************四阶龙格-库塔法**********************/
    double C_R_K_M(double x)
    {
    	if (fabs(x - low) < 1e-6)	//若x==0   第一个初值给定
    		return 1;
    	else {
    		double k1 = C_R_K_M(x - h) - 2 * (x - h) / C_R_K_M(x - h);
    		double k2 = C_R_K_M(x - h) + h / 2.0 * k1 - (2 * (x - h) + h) / (C_R_K_M(x - h) + h / 2.0 * k1);
    		double k3 = C_R_K_M(x - h) + h / 2.0 * k2 - (2 * (x - h) + h) / (C_R_K_M(x - h) + h / 2.0 * k2);
    		double k4 = C_R_K_M(x - h) + h * k3  - (2 * (x - h) + h) / (C_R_K_M(x - h) + h  * k3);
    		return C_R_K_M(x - h) + h / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4);
    	}
    }
    
    
    
    /*************************4、格式化输出函数**********************/
    void prit(double (*p)(double x))      
    {
    	int count=0;
    	double i;
    	for (i = low; high - i >0 ; i += h) {
    		printf("y(%0.3lf)=%0.5lf\t", i, p(i));
    		count++;
    		if (count % 4 == 0)
    			printf("\n");
    
    	}
    	printf("y(%0.3lf)=%0.5lf\t", i, p(i));		//high==i时
    	printf("\n\n\n");
    }
    
    
    int main()
    {
    	printf("解析解格式化输出:\n");
    	prit(AnalyticalSolution);
    	
    	printf("显示欧拉法格式化输出:\n");
    	prit(EulerMethod);
    
    	printf("隐示欧拉法格式化输出:\n");
    	prit(ImplicitEulerMethod_2);
    	
    	printf("两步欧拉法格式化输出:\n");
    	prit(DoubleStepEuler);
    	
    	printf("四阶龙格-库塔法格式化输出:\n");
    	prit(C_R_K_M);
    	return 0;
    }
    

    输出结果:
    在这里插入图片描述

    展开全文
  • 文章目录前言1 常微分方程1.1 常微分方程的概念1.2 常微分方程数值求解的一般概念2 常微分方程数值求解函数3 刚性问题结语 前言 1 常微分方程 1.1 常微分方程的概念 1.2 常微分方程数值求解的一般概念 求解常微分...
  • 专业 序号 姓名 日期 实验 3 常微分方程数值解 实验目的 1掌握用 MATLAB 微分方程初值问题数值解的方法 2通过实例学习微分方程模型解决简化的实际问题 3了解欧拉方法和龙格库塔方法的基本思想 实验内容 用欧拉方法...
  • Euler 法求解常微分方 阿达姆斯预测校正方法求解常微分方程数值解改进的欧拉法求解常微分方程龙格库塔方法求解常微分方程数值解
  • 1.最简模型两点边值问题 2.一般模型二阶常微分方程 3.二维椭圆型方程 4.一维抛物线方程 5.二维非齐次热传导方程
  • 1_基于MATLAB求解常微分方程数值解和解析解.ppt数学建模
  • 耦合常微分方程数值解。 I&EC 基础。 7:3,1968 年。第 514-517 页 测试问题是一个使用 Stefan-Maxwell 方程和两个 Dirichlet 边界条件的三元扩散问题,来自 Ross Taylor & R. Krishna 撰写的“多组分传质”中的第...
  • matlab求解常微分方程数值解

    千次阅读 2017-08-16 19:06:27
    求解如下常微分方程组,其中有s和i两个变量 首先声明这个函数文件 function y=infect(t,x) lamp=11; u=3; y=[lamp*x(1)*x(2)-u*x(1),-lamp*x(1)*x(2)]'; x(1)代表i x(2)代表s 利用ode45对其求数值解 ode...
  • 数值分析实验课代码,包含多种插值方法、多种数值积分方法、多种常微分方程数值求解方法等。
  • 用改进欧拉法一阶常微分方程数值解,计算结果精确
  • MATLAB常微分方程数值解法

    千次阅读 2019-09-25 07:01:27
    科学技术中常常要求解常微分方程的定问题,所谓数值解法就是未知函数在一系列离散点处的近似值。 二、实验原理 三、实验程序 1. 尤拉公式程序 四、实验内容 选一可求解的常微分方程的定问题,分别用...
  • 常微分方程数值求解的命令: 常微分方程的数值解,MATLAB的命令格式为: [t,y]=solver('odefun',tspan,y0,options)其中solver选择ode45等函数名,odefun为根据待解方程或方程组编写的m文件名,tspan为自变量的...
  • Matlab求微分方程数值解

    万次阅读 多人点赞 2020-08-10 17:53:00
    一、Matlab中求微分方程数值解函数 [x,y]=solver('f',ts,x0,options) 1)x代表自变量 2)y代表函数值 3)solver代表求解函数,常用的为ode45->函数图像变化较为平稳,ode15s->函数图像存在突变(一...
  • 龙格库塔 数值分析 ode rkf 常微分方程数值解 fortrain 以及相关论文
  • Russell,常微分方程的边界值问题的数值解,SIAM,1995。 版本 0.0.2 我该如何设置? 主代码存储库是自包含的,只需要 . 运行Pkg.add("BVP")以将包添加到您的 Julia 安装中,并调用Pkg.test("BVP")以运行大量求解器...
  • python初探常微分方程数值解

    千次阅读 2021-01-27 23:36:37
    基于python,采用四阶龙格-库塔算法进行常微分方程数值求解,最终用Apollo卫星运动轨迹作为示例
  • 常微分方程数值解;第五章常微分方程数值解;第五章常微分方程数值解;2数值解的思想;5.2 Euler方法; 18世纪最杰出的数学家之一13岁时入读巴塞尔大学15岁大学毕业16岁获得硕士学位 1727年-1741年20岁-34岁在彼得堡科学...
  • 微分方程中两点边值问题的追赶赴,参考微分方程数值解

空空如也

空空如也

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

求常微分方程的数值解