python-微分方程计算

首先导入数据

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import minimize
import matplotlib.pyplot as pltdata = np.array([[30, 4],[47.2, 6.1],[70.2, 9.8],[77.4, 35.2],[36.3, 59.4],[20.6, 41.7],[18.1, 19],[21.4, 13],[22, 8.3],[25.4, 9.1],[27.1, 7.4],[40.3, 8],[57, 12.3],[76.6, 19.5],[52.3, 45.7],[19.5, 51.1],[11.2, 29.7],[7.6, 15.8],[14.6, 9.7],[16.2, 10.1],[24.7, 8.6]
])

设置参数和定义函数

# Set initial parameters to small random values or default values
initial_guess = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1]# Define ordinary differential equations
def model(y, t, a, b, c, d, e, f):x, z = ydxdt = a * x - b * x * z - e * x**2dydt = -c * z + d * x * z - f * z**2return [dxdt, dydt]#Define error function
def error(params):a, b, c, d, e, f = paramst = np.linspace(0, 1, len(data))  y0 = [data[0, 0], data[0, 1]]  y_pred = odeint(model, y0, t, args=(a, b, c, d, e, f))x_pred, z_pred = y_pred[:, 0], y_pred[:, 1]error_x = np.sum((x_pred - data[:, 0])**2)error_y = np.sum((z_pred - data[:, 1])**2)total_error = error_x + error_yreturn total_error# Use least squares method to fit parameters
result = minimize(error, initial_guess, method='Nelder-Mead')# Get the fitting parameter values
a_fit, b_fit, c_fit, d_fit, e_fit, f_fit = result.x# Simulate the fitted trajectory
t_fit = np.linspace(0, 1, len(data))
y0_fit = [data[0, 0], data[0, 1]]
y_fit = odeint(model, y0_fit, t_fit, args=(a_fit, b_fit, c_fit, d_fit, e_fit, f_fit))
x_fit, z_fit = y_fit[:, 0], y_fit[:, 1]# Plot the fitting results against the original data points
plt.figure(figsize=(10, 6))
plt.scatter(data[:, 0], data[:, 1], label='Raw data', marker='o', color='blue')
plt.plot(x_fit, z_fit, label='Fitting results', linestyle='-', color='red')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()

 

进一步优化

import matplotlib.pyplot as plt# Different optimization methods
methods = ['Nelder-Mead', 'Powell', 'BFGS', 'L-BFGS-B']for method in methods:result = minimize(error, initial_guess, method=method)a_fit, b_fit, c_fit, d_fit, e_fit, f_fit = result.x# Simulate the fitted trajectoryt_fit = np.linspace(0, 1, len(data))y0_fit = [data[0, 0], data[0, 1]]y_fit = odeint(model, y0_fit, t_fit, args=(a_fit, b_fit, c_fit, d_fit, e_fit, f_fit))x_fit, z_fit = y_fit[:, 0], y_fit[:, 1]# Plot the fitting results against the original data pointsplt.figure(figsize=(10, 6))plt.scatter(data[:, 0], data[:, 1], label='Raw data', marker='o', color='blue')plt.plot(x_fit, z_fit, label=f'Fitting results ({method})', linestyle='-', color='red')plt.xlabel('x')plt.ylabel('y')plt.legend()plt.grid()plt.show()print(f"Parameter values fitted using {method} method:a={a_fit}, b={b_fit}, c={c_fit}, d={d_fit}, e={e_fit}, f={f_fit}")

 

Parameter values fitted using Nelder-Mead method:a=1.2173283165346425, b=0.42516102725023064, c=19.726779624261006, d=0.7743814851338301, e=-0.19482192444374966, f=0.37455729849779884

 

Parameter values fitted using Powell method:a=32.49329459442917, b=0.6910719576651617, c=58.98701472032894, d=1.3524516626786816, e=0.47787798383104335, f=-0.5344483269192019

Parameter values fitted using BFGS method:a=1.2171938888848015, b=0.04968374479958104, c=0.9234835772585344, d=0.947268540340848, e=0.010742224447412019, f=0.7985132960108715

 

Parameter values fitted using L-BFGS-B method:a=1.154759061832585, b=0.32168624538800344, c=0.9455699334793284, d=0.9623931795647013, e=0.2936335531513881, f=0.8566315817923148

进一步优化

#Set parameter search range (interval)
bounds = [(-5, 25), (-5, 25), (-5, 25), (-5, 25), (-5, 25), (-5, 25)]# Use interval search to set initial parameter values
result = differential_evolution(error, bounds)a_fit, b_fit, c_fit, d_fit, e_fit, f_fit = result.x# Simulate the fitted trajectory
t_fit = np.linspace(0, 1, len(data))
y0_fit = [data[0, 0], data[0, 1]]
y_fit = odeint(model, y0_fit, t_fit, args=(a_fit, b_fit, c_fit, d_fit, e_fit, f_fit))
x_fit, z_fit = y_fit[:, 0], y_fit[:, 1]# Plot the fitting results against the original data points
plt.figure(figsize=(10, 6))
plt.scatter(data[:, 0], data[:, 1], label='Raw data', marker='o', color='blue')
plt.plot(x_fit, z_fit, label='Fitting results', linestyle='-', color='red')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()print(f"Fitting parameter values:a={a_fit}, b={b_fit}, c={c_fit}, d={d_fit}, e={e_fit}, f={f_fit}")

Fitting parameter values:a=-0.04946907994199101, b=5.5137169943224755, c=0.6909170053541569, d=10.615879287885402, e=-3.1585499451409937, f=18.4110095977882
发现效果竟然变差了
#Set parameter search range (interval)
bounds = [(-0.1, 10), (-0.1, 10), (-0.1, 10), (-0.1,10), (-0.1, 10), (-0.1, 10)]# Use interval search to set initial parameter values
result = differential_evolution(error, bounds)a_fit, b_fit, c_fit, d_fit, e_fit, f_fit = result.x# Simulate the fitted trajectory
t_fit = np.linspace(0, 1, len(data))
y0_fit = [data[0, 0], data[0, 1]]
y_fit = odeint(model, y0_fit, t_fit, args=(a_fit, b_fit, c_fit, d_fit, e_fit, f_fit))
x_fit, z_fit = y_fit[:, 0], y_fit[:, 1]# Plot the fitting results against the original data points
plt.figure(figsize=(10, 6))
plt.scatter(data[:, 0], data[:, 1], label='Raw data', marker='o', color='blue')
plt.plot(x_fit, z_fit, label='Fitting results', linestyle='-', color='red')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()print(f"Fitting parameter values:a={a_fit}, b={b_fit}, c={c_fit}, d={d_fit}, e={e_fit}, f={f_fit}")

 最终优化结果为:

Fitting parameter values:a=10.0, b=0.6320729493793303, c=10.0, d=0.4325244090515547, e=-0.07495645186059174, f=0.18793803443302332

完整代码

创作不易,希望大家多点赞关注评论!!!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/web/25256.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

java线程相关知识点

Java多线程涉及以下几个关键点 1.线程生命周期:理解线程从创建到销毁的各个阶段,包括新建、运行、阻塞、等待、计时等待和终止。 2.线程同步:掌握如何使用synchronized关键字和Lock接口来同步代码,防止数据竞争和死锁。 3.线程间通…

数据分析必备:一步步教你如何用Pandas做数据分析(21)

1、Pandas 可视化 Pandas 可视化是指使用 Pandas 库中的函数和方法来创建数据可视化图表。Pandas 提供了一些基本的绘图功能,例如折线图、柱状图、饼图等,可以通过调用相应的函数来创建这些图表。 2、基本绘图:绘图 Series和DataFrame上的…

预期值与实际值对比

编辑实际值和预期值变量 因为在单独的代码当中,我们先定义了变量str,所以在matcher时传入str参数,但当我们要把这串代码写在testrun当中,改下传入的参数,与excel表做连接 匹配的结果是excel表中的expect结果&#xf…

等级保护政策法规解读:构建网络安全的法律基石

等级保护政策法规解读:构建网络安全的法律基石 引言 等级保护制度作为中国网络安全管理的基石,其政策法规构成了网络运营者履行安全保护义务的法律框架。随着技术的发展和网络安全形势的变化,等级保护政策法规也在不断更新和完善。本文旨在解…

Python的列表和元组之间的区别是?在 Python 中,如何使用列表和元组进行高效的数据操作?

Python 中的列表(List)和元组(Tuple)是两种不同的数据结构,它们有以下主要区别: 可变性: 列表是可变的(Mutable),这意味着你可以在创建列表后添加、删除或更改…

有序二叉树java实现

类实现: package 树;import java.util.LinkedList; import java.util.Queue;public class BinaryTree {public TreeNode root;//插入public void insert(int value){//插入成功之后要return结束方法TreeNode node new TreeNode(value);//如果root为空的话插入if(r…

RK3288 android7.1 实现ota升级时清除用户数据

一,OTA简介(整包,差分包) OTA全称为Over-The-Air technology(空中下载技术),通过移动通信的接口实现对软件进行远程管理。 1. 用途: OTA两种类型最大的区别莫过于他们的”出发点“(我们对两种不同升级包的创建&…

SolidityFoundry 安全审计测试 Delegatecall漏洞

名称:Delegatecall漏洞 描述: 代理合约所有者操纵漏洞,是智能合约设计中的一个缺陷,允许攻击者操纵代理合约所有者。该漏洞允许攻击者操纵代理合约的所有者(这里我们把所有者硬编码为 0xdeadbeef)。漏洞产…

牛客多校Ancestor(lca,集合的lca)

题目描述 NIO is playing a game about trees. The game has two trees A,BA, BA,B each with NNN vertices. The vertices in each tree are numbered from 111 to NNN and the iii-th vertex has the weight viv_ivi​. The root of each tree is vertex 1. Given KKK key n…

PHP实名认证接口开发示例、银行卡实名认证API

在互联网技术多元化、高速的发展下,催生出在挑战中不断奋勇前进的互联网企业。但不能忽视的是,互联网技术的快速迭代也会使部分企业在冲击中败下阵来,面临淘汰的危机。随着O2O、共享经济等新兴商业形式的兴起,企业对实名认证业务的…

如何使用Python中的列表解析(list comprehension)进行高效列表操作

Python中的列表解析(list comprehension)是一种创建列表的简洁方法,它可以在单行代码中执行复杂的循环和条件逻辑。列表解析提供了一种快速且易于阅读的方式来生成新的列表。 以下是一些使用列表解析进行高效列表操作的示例: 1.…

用Python编写自动发送每日电子邮件报告的脚本

为了用 Python 编写自动发送每日电子邮件报告的脚本,你可以使用 smtplib 库来发送电子邮件,使用 email 库来创建电子邮件内容。此外,你可以使用 schedule 库来安排每天发送邮件的任务。以下是一个示例脚本以及如何设置和运行它的指导。 步骤…

JSON如何处理包含特殊字符的字段

在JSON中处理包含特殊字符的字段时,你通常不需要直接处理这些特殊字符,因为JSON格式本身就会对特殊字符进行转义。当你使用编程语言或工具来生成或解析JSON时,这些转义通常是自动处理的。 然而,如果你需要手动处理或理解这些转义…

华为策略流控

以下脚本仅做参考,具体IP地址和接口请按照现场实际情况写入。 [Huawei]acl 3001 [Huawei-acl-adv-3001]rule permit ip source 192.168.1.10 0.0.0.0 destination 192.168.2.10 0.0.0.0 //匹配需要做测试的源和目标地址 [Huawei-acl-adv-3001]rule permit ip sour…

[AIGC] CompletableFuture的重要方法有哪些?

CompletableFuture具有多种方法&#xff0c;使其成为异步编程的强大工具。在这里&#xff0c;我们将介绍一些最重要和常用的方法&#xff1a; CompletableFuture<T> supplyAsync(Supplier<T> supplier): 使用ForkJoinPool.commonPool()作为线程池来异步执行Suppile…

力扣2781.最长合法子字符串的长度

力扣2781.最长合法子字符串的长度 将字符串数组存入哈希表 枚举所有右端点反向遍历子串在哈希表中找所有以i为右端点的字符串若找到相同子串 更新j k 1 class Solution {public:int longestValidSubstring(string word, vector<string>& forbidden) {unordered_…

【马琴绿绮】马维衡古琴之马氏汉风 明代杉木制;周身髹朱红色漆

【马琴绿绮式】马维衡古琴之马氏汉风 明代杉木制&#xff1b;琴体周身髹朱红色漆&#xff0c;鹿角霜灰胎&#xff1b;形体壮硕、风格高古&#xff1b;音色松透、浑厚&#xff0c;音质纯净&#xff0c;按弹舒适&#xff0c;手感丝滑。

C++ 课堂实验 读取a.txt中文本,统计文本中字母数量

题目描述:读取a.txt中文本&#xff0c;统计文本中字母数量。 相关知识&#xff08;略&#xff09; 编程要求 根据提示&#xff0c;在右侧编辑器Begin-End处补充代码&#xff0c;完成本关要求。 测试说明 输入 读取a.txt读入文本 如&#xff1a; abc abc 输出 输出文本中字母数…

Effective Java 2 遇到多个构造器参数时要考虑使用构建器

第2个经验法则&#xff1a;用遇到多个构造器参数时要考虑使用构建器&#xff08;consider a builder when faced with many constructor parameters&#xff09; 上一条讨论了静态工厂相对于构造器来说有五大优势。但静态工厂和构造器有个共同的局限性:它 们都不能很好地扩展到…

华为坤灵路由器初始化的几个坑,含NAT配置

1、aaa密码复杂度修改&#xff1a; #使能设备对密码进行四选三复杂度检查功能。 <HUAWEI>system-view [HUAWEI]aaa [HUAWEI-aaa]local-aaa-user password policy administrator [HUAWEI-aaa-lupp-admin]password complexity three-of-kinds 2、本地用户名长度必须大…