当前位置: 首页 > news >正文

Fortran 微分方程求解 --ODEPACK

最近涉及到使用Fortran对微分方程求解,我们知道MATLAB已有内置的函数,比如ode家族,ode15s,对应着不同的求解办法。通过查看odepack的官方文档,我尝试使用了dlsode求解刚性和非刚性常微分方程组。

首先是github网址:https://github.com/jacobwilliams/odepack

具体使用办法:

1.我使用的是vs2022,比较简单的用法就是把,src文件夹所有的文件复制到和项目一个文件夹即可,将M_odepack.f90文件放入到项目中,这样就可以用了。

2.在使用前要use m_odepack

3.这里以官方文档中的例子为例:

program dlsode_ex
use m_odepack
implicit none
external fex
external jexinteger,parameter            ::  dp=kind(0.0d0)
real(kind=dp),dimension(3)   ::  atol,y
integer                      ::  iopt,iout,istate,itask,itol,liw,lrw,mf,neq
integer,dimension(23)        ::  iwork
real(kind=dp)                ::  rtol,t,tout
real(kind=dp),dimension(58)  ::  rworkneq = 3y(1) = 1.D0y(2) = 0.D0y(3) = 0.D0t = 0.D0tout = .4D0itol = 2rtol = 1.D-4atol(1) = 1.D-6atol(2) = 1.D-10atol(3) = 1.D-6itask = 1istate = 1iopt = 0lrw = 58liw = 23mf = 21do iout = 1,12call dlsode(fex,[neq],y,t,tout,itol,[rtol],atol,itask,istate,iopt,   && rwork,lrw,iwork,liw,jex,mf)write (6,99010) t,y(1),y(2),y(3)99010 format (' At t =',d12.4,'   y =',3D14.6)if ( istate<0 ) thenwrite (6,99020) istate99020 format (///' Error halt.. ISTATE =',i3)stop 1elsetout = tout*10.D0endifenddowrite (6,99030) iwork(11),iwork(12),iwork(13)99030 format (/' No. steps =',i4,',  No. f-s =',i4,',  No. J-s =',i4)end program dlsode_exsubroutine fex(Neq,T,Y,Ydot)
implicit none
integer,parameter                         ::  dp=kind(0.0d0)integer                                   ::  Neq
real(kind=dp)                             ::  T
real(kind=dp),intent(in),dimension(3)     ::  Y
real(kind=dp),intent(inout),dimension(3)  ::  YdotYdot(1) = -.04D0*Y(1) + 1.D4*Y(2)*Y(3)Ydot(3) = 3.D7*Y(2)*Y(2)Ydot(2) = -Ydot(1) - Ydot(3)
end subroutine fexsubroutine jex(Neq,T,Y,Ml,Mu,Pd,Nrpd)
implicit noneinteger,parameter                              ::  dp=kind(0.0d0)
integer                                        ::  Neq
real(kind=dp)                                  ::  T
real(kind=dp),intent(in),dimension(3)          ::  Y
integer                                        ::  Ml
integer                                        ::  Mu
real(kind=dp),intent(inout),dimension(Nrpd,3)  ::  Pd
integer,intent(in)                             ::  NrpdPd(1,1) = -.04D0Pd(1,2) = 1.D4*Y(3)Pd(1,3) = 1.D4*Y(2)Pd(2,1) = .04D0Pd(2,3) = -Pd(1,3)Pd(3,2) = 6.D7*Y(2)Pd(2,2) = -Pd(1,2) - Pd(3,2)
end subroutine jex

一些变量意义具体看文档说明:https://jacobwilliams.github.io/odepack/proc/dlsode.html

其中,假设n是方程个数,

y:是初值,数组,y(n)

atol:每个方程的绝对误差,数组,atol(n)

t:输入的初始点,tout是下一个点。

mf:是求解方法,其中如果等于21,24需要使用者自己提供雅各比矩阵,如示例代码中jex函数中那样,如果等于10,22,25则不需要自己写,但是jex函数还是需要定义,就是函数框架,函数名,变量声明就可。

fex函数:写的就是你的微分方程组

另外,

 rwork,iwork也是两个一维数组,大小如图所示。

以及,

lrw = 22 +  9*NEQ + NEQ**2
liw = 20 + NEQ

整体使用的逻辑就是先设置t值,然后设置循环,tout不断累加,下次循环就使用上次计算得到的新y值以及tout进行迭代计算。

istate是用于输入和输出以指定计算状态的索引,要注意的是如果istate选择2,或3需要在第一次循环中等于1,初始化,到了第二次循环开始才赋值为2或3。

相关文章:

Fortran 微分方程求解 --ODEPACK

最近涉及到使用Fortran对微分方程求解&#xff0c;我们知道MATLAB已有内置的函数&#xff0c;比如ode家族&#xff0c;ode15s&#xff0c;对应着不同的求解办法。通过查看odepack的官方文档&#xff0c;我尝试使用了dlsode求解刚性和非刚性常微分方程组。 首先是github网址&am…...

8路光栅尺磁栅尺编码器或16路高速DI脉冲信号转Modbus TCP网络模块 YL99-RJ45

特点&#xff1a; ● 光栅尺磁栅尺解码转换成标准Modbus TCP协议 ● 高速光栅尺磁栅尺4倍频计数&#xff0c;频率可达5MHz ● 模块可以输出5V的电源给光栅尺或传感器供电 ● 支持8个光栅尺同时计数&#xff0c;可识别正反转 ● 可以设置作为16路独立DI高速计数器 ● 可网…...

【Python】函数

None类型 思考&#xff1a;若函数没有使用return语句返回数据&#xff0c;那么函数有返回值吗&#xff1f; 答&#xff1a;实际上是有的&#xff0c;Python中有一个特殊的字面量None&#xff0c;其类型是<class ‘NoneType’>&#xff0c;无返回值的函数&#xff0c;实…...

centos安装MySQL 解压版完整教程(按步骤傻瓜式安装

一、卸载系统自带的 Mariadb 查看&#xff1a; rpm -qa|grep mariadb 卸载&#xff1a; rpm -e --nodeps mariadb-libs-5.5.68-1.el7.x86_64 二、卸载 etc 目录下的 my.cnf 文件 rm -rf /etc/my.cnf 三、检查MySQL是否存在 有则先删除 #卸载mysql服务以及删除所有mysql目录 #没…...

【后端速成 Vue】第一个 Vue 程序

1、为什么要学习 Vue&#xff1f; 为什么使用 Vue? 回想之前&#xff0c;前后端交互的时候&#xff0c;前端收到后端响应的数据&#xff0c;接着将数据渲染到页面上&#xff0c;之前使用的是 JavaScript 或者 基于 JavaScript 的 Jquery&#xff0c;但是这两个用起来还是不太…...

Macbook pro M1 安装Ubuntu教程

先讲下心路历程 由于版主最近刚切换到Mac&#xff0c;所以在安装的时候一上手就选择了virutalbox&#xff0c;结果报错“The installer has detected an unsupported architecture. VirtualBox only runs on the amd64 architecture.” 后来去Reddit论坛上一看&#xff0c;才知…...

前端console.log打印内容与后端请求返回数据不一致

后端传值num0 前端打印num1 ,如图&#xff0c;console.log后台显示的数据与展开后不一致 造成该问题原因是深拷贝与浅拷贝的问题。 var obj JSON.parse(JSON.stringify(res)) 修改后打印 正常...

SQL入门:多表查询

SQL&#xff0c;或者说结构化查询语言(Structured Query Language)&#xff0c;是用于管理和操作关系型数据库的标准语言。在本篇文章中&#xff0c;我们将重点介绍SQL中的多表查询&#xff0c;这是一种强大的工具&#xff0c;可以帮助我们从多个相关的表格中获取数据。 数据库…...

【C++】进一步认识模板

&#x1f3d6;️作者&#xff1a;malloc不出对象 ⛺专栏&#xff1a;C的学习之路 &#x1f466;个人简介&#xff1a;一名双非本科院校大二在读的科班编程菜鸟&#xff0c;努力编程只为赶上各位大佬的步伐&#x1f648;&#x1f648; 目录 前言一、非类型模板参数二、模板的特…...

Mysql Oracle 区别

1. oracle select *, id需要在星号前加别名&#xff0c;mysql则不需要 mysql语法&#xff1a; select *, id from xin_student_t;oracle语法&#xff1a; select st.*, st.id from xin_student_t st;2. oracle表定义了别名&#xff0c;在查询时可以不用别名指定字段&#xf…...

华为OD-第K长的连续字母字符串长度

题目描述 给定一个字符串&#xff0c;只包含大写字母&#xff0c;求在包含同一字母的子串中&#xff0c;长度第 k 长的子串的长度&#xff0c;相同字母只取最长的那个子串。 代码实现 # coding:utf-8 # 第K长的连续字母字符串长度 # https://www.nowcoder.com/discuss/353150…...

【编程题】有效三角形的个数

文章目录 一、题目二、算法讲解三、题目链接四、补充 一、题目 给定一个包含非负整数的数组 nums &#xff0c;返回其中可以组成三角形三条边的三元组个数。 示例1&#xff1a; 输入: nums [2,2,3,4] 输出: 3 **解释:**有效的组合是: 2,3,4 (使用第一个 2) 2,3,4 (使用第二个 …...

【mysql是怎样运行的】-EXPLAIN详解

文章目录 1.基本语法2. EXPLAIN各列作用1. table2. id3. select_type4. partitions5. type 1.基本语法 EXPLAIN SELECT select_options #或者 DESCRIBE SELECT select_optionsEXPLAIN 语句输出的各个列的作用如下&#xff1a; 列名描述id在一个大的查询语句中每个SELECT关键…...

数据结构例题代码及其讲解-链表

链表 单链表的结构体定义及其初始化。 typedef struct LNode {int data;struct LNode* next; }LNode, *LinkList;①强调结点 LNode *p; ②强调链表 LinkList p; //初始化 LNode* initList() {//定义头结点LNode* L (LNode*)malloc(sizeof(LNode));L->next NULL;return …...

[Open-source tool] 可搭配PHP和SQL的表單開源工具_Form tools(1):簡介和建置

Form tools是一套可搭配PHP和SQL的表單開源工具&#xff0c;可讓開發者靈活運用&#xff0c;同時其有數個表單模板和應用模組供挑選&#xff0c;方便且彈性。Form tools已開發超過20年&#xff0c;為不同領域的需求者或開發者提供一個自由和開放的平台&#xff0c;使他們可建構…...

移动数据业务价值链的整合

3G 时代移动数据业务开发体系的建立和发展&#xff0c;要求运营商从封闭、统一的业 务形态、单一提供业务&#xff0c;向开放的、个性化多元化的业务体系以及多方合作参与提 供业务的方向发展&#xff0c;不可避免的使通信价值链不断延长和升级&#xff0c;内容提供商、服务 …...

合并两个链表

题目描述 将两个升序链表合并为一个新的 升序 链表并返回。新链表是通过拼接给定的两个链表的所有节点组成的。 比如以下例子&#xff1a; 题目接口&#xff1a; /*** Definition for singly-linked list.* struct ListNode {* int val;* ListNode *next;* ListN…...

测试框架pytest教程(9)跳过测试skip和xfail

skip无条件跳过 使用装饰器 pytest.mark.skip(reason"no way of currently testing this") def test_example(faker):print("nihao")print(faker.words()) 方法内部调用 满足条件时跳过 def test_example():a1if a>0:pytest.skip("unsupported …...

HTML <textarea> 标签

实例 <textarea rows="3" cols="20"> 收拾收拾 </textarea>定义和用法 <textarea> 标签定义多行的文本输入控件。 文本区中可容纳无限数量的文本,其中的文本的默认字体是等宽字体(通常是 Courier)。 可以通过 cols 和 rows 属性来…...

探索图结构:从基础到算法应用

文章目录 理解图的基本概念学习图的遍历算法学习最短路径算法案例分析&#xff1a;使用 Dijkstra 算法找出最短路径结论 &#x1f389;欢迎来到数据结构学习专栏~探索图结构&#xff1a;从基础到算法应用 ☆* o(≧▽≦)o *☆嗨~我是IT陈寒&#x1f379;✨博客主页&#xff1a;I…...

大模型赋能金融底稿搜索:告别大海捞针,实现高效精准合规管理!

文章主要介绍了达观数据利用大模型技术升级其底稿搜索产品&#xff0c;为金融行业带来革命性的变化。传统底稿搜索存在关键词匹配局限、非结构化文件解析困难、溯源关联不便和合规风险高等问题。达观数据通过深度语义理解、全格式解析兼容、智能要素抽取、全链路溯源关联和开箱…...

告别向日葵和TeamViewer!用你家路由器自带的DDNS功能,免费搭建Windows远程桌面(保姆级教程)

告别第三方远程工具&#xff1a;用路由器DDNS解锁Windows远程桌面全速体验 每次打开向日葵或TeamViewer时&#xff0c;那个转圈加载的进度条是否让你眉头紧锁&#xff1f;当免费版突然弹出"会话时长已达上限"的提示时&#xff0c;是否恨不得砸键盘&#xff1f;作为常…...

51单片机项目避坑:用ADC0804读PT100信号,你的滤波和标度变换做对了吗?(附源码分析)

51单片机PT100温度检测实战&#xff1a;从ADC采样到标度变换的完整设计解析 在工业温度测量领域&#xff0c;PT100凭借其优异的线性度和稳定性成为首选传感器之一。不同于常见的DS18B20数字温度传感器&#xff0c;PT100需要配合精密信号调理电路和AD转换器才能实现准确测量。本…...

自动驾驶、无人机导航都离不开它:卡尔曼滤波在机器人SLAM中的实战调参心得

自动驾驶与无人机导航中的卡尔曼滤波实战&#xff1a;SLAM系统调参进阶指南 卡尔曼滤波算法自1960年问世以来&#xff0c;已成为机器人定位与导航领域不可或缺的核心技术。无论是自动驾驶汽车的精准定位&#xff0c;还是无人机在复杂环境中的自主飞行&#xff0c;亦或是工业机器…...

基于主从博弈的主动配电网阻塞管理探索

基于主从博弈的主动配电网阻塞管理 首先&#xff0c;在日前市场中&#xff0c;LA&#xff08;负荷聚合商&#xff09;根据历史数据预测次日向上级电网购电的电价信息和预测分布式电源(燃气轮机)出力、风电场出力信息&#xff0c;同时考虑事前与用户签订协议的可中断负荷&#x…...

ArXiv:为何大模型无法拥有意识|Erik Hoel

导语当AI能流畅谈论“自我感受”&#xff0c;当Anthropic赋予Claude“对话退出权”&#xff0c;我们是否可以说它有意识&#xff1f;2026年初&#xff0c;神经科学家Erik Hoel在ArXiv发布论文《大语言模型意识证伪&#xff1a;持续学习对意识存在的必要性》&#xff08;A Dispr…...

Redis未授权访问漏洞实战:从SSH公钥到反弹shell的5种利用方式详解

Redis未授权访问漏洞深度攻防&#xff1a;5种高阶利用与防御方案 Redis作为高性能键值数据库&#xff0c;其未授权访问漏洞长期位居企业安全风险Top 10。本文将突破常规教程框架&#xff0c;从攻击者视角剖析5种实战利用手法&#xff0c;同时提供企业级防御方案。不同于基础复现…...

Vue中实现动态标签页的切换优化与状态管理

1. 动态标签页的核心需求与实现思路 在后台管理系统这类多页面应用中&#xff0c;动态标签页几乎是标配功能。想象一下你正在使用某电商后台&#xff0c;同时开着商品管理、订单处理和用户分析三个页面&#xff0c;这时候标签页的流畅切换和状态保持就显得尤为重要。 我经历过一…...

在 React 中,useRef、ref 属性以及 forwardRef 是处理“引用”(访问 DOM 节点或组件实例)的核心概念

在 React 中&#xff0c;useRef、ref 属性以及 forwardRef 是处理“引用”&#xff08;访问 DOM 节点或组件实例&#xff09;的核心概念。它们经常一起使用&#xff0c;但职责完全不同。以下是它们的核心区别、使用方法及组合示例&#xff1a;1. 核心概念与区别特性ref (属性)u…...

飞行器设计避坑指南:盘点那些影响气动效率的‘隐形杀手’(从摩擦阻力到干扰阻力)

飞行器设计避坑指南&#xff1a;盘点那些影响气动效率的‘隐形杀手’ 记得第一次参加大学生飞行器设计竞赛时&#xff0c;我们的团队花了整整三个月打造了一架翼展两米的固定翼无人机。试飞当天&#xff0c;看着它摇摇晃晃地起飞&#xff0c;却在爬升阶段突然失速坠毁&#xff…...