当前位置: 首页 > 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…...

3分钟掌握Umi-OCR插件:打造你的专属文字识别工具箱

3分钟掌握Umi-OCR插件&#xff1a;打造你的专属文字识别工具箱 【免费下载链接】Umi-OCR_plugins Umi-OCR 插件库 项目地址: https://gitcode.com/gh_mirrors/um/Umi-OCR_plugins 还在为不同场景下的文字识别需求而烦恼吗&#xff1f;Umi-OCR插件库为你提供了完美的解决…...

智能家居控制中心:OpenClaw桥接Qwen3-32B-Chat与HomeAssistant

智能家居控制中心&#xff1a;OpenClaw桥接Qwen3-32B-Chat与HomeAssistant 1. 为什么需要AI驱动的家居控制中心 去年冬天的一个深夜&#xff0c;我被空调异常制热的噪音惊醒。摸黑在手机APP上反复调整参数无果后&#xff0c;突然意识到&#xff1a;如果有个能理解自然语言的智…...

HTML网页元素中的图片和超链接

哈哈哈&#xff0c;又来更新我这一周里面新学的web前端开发技术啦&#xff01;今天我将与大家分享网页元素中的图片和超链接。一.图像的应用HTML中加入图片有3种不同的路径&#xff1a;1.绝对路径&#xff1a;是指互联网上唯一且完整的地址&#xff0c;用来精准定位资源。绝对路…...

基于Matlab的11种图像清晰度评价指标:直接可运行,联系我

基于matlab图像清晰度评价指标。 一共11种。 程序已调通&#xff0c;可直接运行。 需要直接联系。 基于matlab图像清晰度评价指标。 一共11种。 程序已调通&#xff0c;可直接运行。 需要直接联系。 图像剃度的清晰度评价(EOG, Roberts, Tenengrad, Brenner,Variance, Laplace,…...

隐私保护方案:OpenClaw+GLM-4.7-Flash本地化处理敏感数据

隐私保护方案&#xff1a;OpenClawGLM-4.7-Flash本地化处理敏感数据 1. 为什么需要本地化处理敏感数据&#xff1f; 去年我帮一位做财务咨询的朋友处理季度报表时&#xff0c;遇到了一个棘手问题。他需要分析上百份包含客户银行流水、身份证号等信息的Excel文件&#xff0c;但…...

避坑指南:关系数据库设计中90%人会犯的完整性约束错误(附真实案例)

避坑指南&#xff1a;关系数据库设计中90%人会犯的完整性约束错误&#xff08;附真实案例&#xff09; 在电商大促期间&#xff0c;某平台突然出现大量"幽灵订单"——用户支付成功后订单消失&#xff0c;而库存却异常扣减。技术团队紧急排查发现&#xff0c;问题根源…...

nli-distilroberta-base环境配置:Ubuntu/CentOS下Python依赖与端口映射设置

nli-distilroberta-base环境配置&#xff1a;Ubuntu/CentOS下Python依赖与端口映射设置 1. 项目概述 nli-distilroberta-base是一个基于DistilRoBERTa模型的自然语言推理(NLI)Web服务&#xff0c;专门用于判断两个句子之间的逻辑关系。这个轻量级模型保留了RoBERTa-base模型9…...

零成本实现3D模型跨平台迁移:Blender到Unreal Engine的无缝解决方案

零成本实现3D模型跨平台迁移&#xff1a;Blender到Unreal Engine的无缝解决方案 【免费下载链接】bl_datasmith Blender addon to export UE4 Datasmith format 项目地址: https://gitcode.com/gh_mirrors/bl/bl_datasmith 你是否曾遇到这样的困境&#xff1a;在Blender…...

ROS2 MoveIt2实战:如何让虚拟机械臂‘看懂’并抓取YOLOv8 OBB识别的物体?

ROS2 MoveIt2与YOLOv8 OBB深度集成&#xff1a;构建高精度虚拟抓取系统的核心技术解析 当机械臂遇上计算机视觉&#xff0c;一场关于精准控制的交响乐就此展开。本文将带您深入探索如何利用YOLOv8 OBB&#xff08;Oriented Bounding Box&#xff09;的朝向感知能力&#xff0c;…...

Vue-Sonner:面向现代Vue应用的高性能Toast通知架构解析

Vue-Sonner&#xff1a;面向现代Vue应用的高性能Toast通知架构解析 【免费下载链接】vue-sonner &#x1f514; An opinionated toast component for Vue. 项目地址: https://gitcode.com/gh_mirrors/vu/vue-sonner 在当今快节奏的Web应用开发中&#xff0c;实时反馈机制…...