解双曲型非线性方程的Harden-Yee算法(TVD格式)

解双曲型非线性方程的Harden-Yee算法
先贴代码,教程后面有空再写

import matplotlib
import math
matplotlib.use('TkAgg')
import numpy as np  
import matplotlib.pyplot as plt
def Phiy(yy,epsi):#phi(y)if abs(yy) >= epsi:phiyy = abs(yy)else:phiyy = (yy*yy + epsi*epsi)/2.*epsireturn phiyy
def Deltau(u2, u1):#\delta u(n,i+1/2),u2表示右边的return u2 - u1def Gi_0(u2,u1,u0):result = Minmod((u1-u0),(u2-u1))return resultdef Minmod(aa,bb):if aa*bb<=0.:result = 0.else:result = np.sign(aa)*min(abs(aa),abs(bb))return resultdef Hartenyee(u,x,t):corrant = (t[2] - t[1])/( x[2]-x[1])#corrant 数epsi=0.3for j in range(0, t.size-1):# j 表示t方向,t从零开始u[j+1,1]=Lax_Wendroff(u[j,0],u[j,1],u[j,2],corrant)#赋值x1print(j)for i in range(2, x.size-2): # i 表示x 方向if Deltau(u[j,i+1],u[j,i]) == 0.:alphai_1plus2 = u[j,i]beta_1plus2 = 0.else:alphai_1plus2 = (u[j,i+1] + u[j,i])/2.#E2-E1beta_1plus2 = (Gi_0(u[j,i+2],u[j,i+1],u[j,i]) - Gi_0(u[j,i+1],u[j,i],u[j,i-1]))/(u[j,i+1]-u[j,i])if Deltau(u[j,i],u[j,i-1]) == 0.:alphai_1minus2 = u[j,i]beta_1minus2 = 0.else:alphai_1minus2 = (u[j,i] + u[j,i-1])/2.#E2-E1beta_1minus2 = (Gi_0(u[j,i+1],u[j,i],u[j,i-1]) - Gi_0(u[j,i],u[j,i-1],u[j,i-2]))/(u[j,i]-u[j,i-1])phini_plus12 = Gi_0(u[j,i+2],u[j,i+1],u[j,i])+Gi_0(u[j,i+1],u[j,i],u[j,i-1])-Phiy(alphai_1plus2+beta_1plus2,epsi)*(u[j,i+1]-u[j,i])phini_minus12 = Gi_0(u[j,i+1],u[j,i],u[j,i-1])+Gi_0(u[j,i],u[j,i-1],u[j,i-2])-Phiy(alphai_1minus2+beta_1minus2,epsi)*(u[j,i]-u[j,i-1])u[j+1,i] = u[j,i] - corrant/2.*(u[j,i+1]*u[j,i+1]/2.-u[j,i-1]*u[j,i-1]/2.)-corrant/2*(phini_plus12-phini_minus12)return u
def Lax_Wendroff(u0,u1,u2,corrant):   #lax——Wendroff格式dE1 = u2*u2/4.-u0*u0/4.dE2 = u2*u2/4.-u1*u1/4.dE3 = u1*u1/4.-u0*u0/4.result = u1 - corrant/2.*dE1 + corrant*corrant/2.*((u1+u2)/2.*dE2-(u1+u0)/2.*dE3) return resultdef Plot(x, t, result,title):plt.figure()plt.plot(x, result[0,:])y = [t[int(t.size/5)], t[int(t.size/4)], t[int(t.size/3)], t[int(t.size/2)],t[int(t.size/1.1)]]labels = ['t='f'{num:.2f}' for num in y]#将变量转换为字符串plt.plot(x, result[int(t.size/5),:], label=labels[0])plt.plot(x, result[int(t.size/4),:],label=labels[1])plt.plot(x, result[int(t.size/3),:],label=labels[2])plt.plot(x, result[int(t.size/2),:],label=labels[3])plt.plot(x, result[int(t.size/1.1),:],label=labels[4])plt.legend()plt.xlabel('x')plt.ylabel('t')plt.title(title)plt.show()plt.close()plt.figure()plt.contourf(x, t, result, 50, cmap = 'jet')plt.colorbar()plt.savefig('CN.png')plt.xlabel('x')plt.ylabel('t')plt.title(title)plt.show()plt.close()return 0x = np.linspace(0,4,400)
t = np.linspace(0,4.0,400)
u = np.zeros((t.size,x.size),dtype=float)#注意,这里u不是必须二维,我用二维主要为了测试和画图方便,当t比较大的时候会占用大量内存
u[0,0:int(x.size/2)] = 1.0
u[:,0] = 1.0
u[:,-1] = 0.0
corrant = (t[2] - t[1])/( x[2]-x[1])#corrant 数
'''#Lax_Wendroff 方法做对比
for j in range (1, t.size-1):print(j)for i in range (1,x.size-1):u[j,i] = Lax_Wendroff(u[j-1,i-1],u[j-1,i],u[j-1,i+1],corrant)
'''u2=Hartenyee(u,x,t)
Plot(x,t,u2,'Harten-Yee solver')

在这里插入图片描述
在这里插入图片描述

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

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

相关文章

国外在线教育系统源码,知识付费课程录制流程是什么样?

无论是从信息内容优化&#xff0c;还是知识产权保护的角度来看&#xff0c;“内容付费”都是一个有力的抓手&#xff0c;并且也一定是未来互联网发展的一个重要方向。与此同时&#xff0c;互联网技术的不断进步降低了“内容付费”的使用门槛&#xff0c;越来越多的人企图搭上这…

2024年华为OD机试真题-螺旋数字矩阵-Java-OD统一考试(C卷D卷)

题目描述: 疫情期间,小明隔离在家,百无聊赖,在纸上写数字玩。他发明了一种写法: 给出数字个数n和行数m(0 < n ≤ 999,0 < m ≤ 999),从左上角的1开始,按照顺时针螺旋向内写方式,依次写出2,3...n,最终形成一个m行矩阵。 小明对这个矩阵有些要求: 1.每行数字的…

2024年记录

kong docker 安装: kong docker 安装

C++ VScode: launch: program ...... dose not exist

VScode: launch: program … dose not exist 介绍 参考VS Code 配置 C/C 编程运行环境&#xff08;保姆级教程&#xff09;教程配置了VSCode。在配置launch.json适用多个.c 文件编译时&#xff0c;弹出下面错误。 原因和解决方法 是task.json 默认配置的问题。 默认的 cwd参…

【深度学习Labelme】使用Segment Anything Model (SAM)快速打标,labelme多边形转yolo txt框看看对不对

文章目录 windows安装环境打开labelme自动保存勾选上&#xff0c;保存图片数据不要勾选选SAM精准模型&#xff0c;然后打开图片路径&#xff0c;然后点击创建AI多边形&#xff1a;鼠标点击确认物体控制点&#xff0c;确认完成后&#xff0c;双击鼠标完成选取&#xff0c;并给上…

利用OpenShift的ImageStream部署临时版本

公司是港企&#xff0c;项目都部署在OpenShift上统一管理&#xff0c;因为运行环境为香港网络(外网)&#xff0c;配置、中间件等大陆无法直接访问联通。因此在大陆开发时&#xff0c;测试是个很大的问题。为了避免往Git上频繁提交未确定可用的版本&#xff0c;选择用利用OpenSh…

知识付费系统方面的源码,如何提高培训成果转化的可能性?重心应放在哪?

培训成果转化对于一个教育机构来说是非常关键的&#xff0c;因为不管是学员的成果还是员工的培训&#xff0c;它的目的性都非常的关键&#xff0c;但是如何提高成果转化的可能性?这个在很大程度上来说&#xff0c;可能大家还存有疑问。 一、提高培训成果转化的可能性 成人学习…

嵌入式学习70-复习(wireshark使用和http协议)

--------------------------------------------------------------------------------------------------------------------------------- wireshark 1.sudo wireshark 2.选择 any &#xff0c; 3.搜索 http/tcp 54 为 发送的数据包 58 回复的数据包 请求报文 请求报文…

【NLP练习】使用seq2seq实现文本翻译

使用seq2seq实现文本翻译 &#x1f368; 本文为&#x1f517;365天深度学习训练营 中的学习记录博客&#x1f356; 原作者&#xff1a;K同学啊 from __future__ import unicode_literals, print_function, division from io import open import unicodedata import string impo…

05、Kafka 操作命令

05、Kafka 操作命令 1、主题命令 &#xff08;1&#xff09;创建主题 kafka-topics.sh --create --bootstrap-server 192.168.135.132:9092,192.168.135.133:9092,192.168.135.134:9092 --topic test1 --partitions 4 --replication-factor 3–bootstrap-server&#xff1a;…

Gradient发布支持100万token的Lllama3,上下文长度从8K扩展到1048K

前言 近日Gradient公司在Crusoe Energy公司的算力支持下&#xff0c;开发了一款基于Llama-3的大型语言模型。这款新模型在原Llama-3 8B的基础上&#xff0c;将上下文长度从8000 token大幅扩展到超过104万token。 这一创新性突破&#xff0c;展现了当前SOTA大语言模型在长上下…

类和对象(上篇)

面向对象和面向过程 C语言是面向过程的&#xff0c;关注的是过程&#xff0c;分析出求解问题的步骤&#xff0c;通过函数调用逐步解决问题。 C是基于面向对象的&#xff0c;关注的是对象&#xff0c;将一件事情拆分成不同的对象&#xff0c;靠对象之间的交互完成。 类的引入 在…

【17-Ⅰ】Head First Java 学习笔记

HeadFirst Java 本人有C语言基础&#xff0c;通过阅读Java廖雪峰网站&#xff0c;简单速成了java&#xff0c;但对其中一些入门概念有所疏漏&#xff0c;阅读本书以弥补。 第一章 Java入门 第二章 面向对象 第三章 变量 第四章 方法操作实例变量 第五章 程序实战 第六章 Java…

024.反转链表

给定单链表的头节点 head &#xff0c;请反转链表&#xff0c;并返回反转后的链表的头节点。 示例 1&#xff1a; 输入&#xff1a;head [1,2,3,4,5] 输出&#xff1a;[5,4,3,2,1]示例 2&#xff1a; 输入&#xff1a;head [1,2] 输出&#xff1a;[2,1]示例 3&#xff1a; 输…

《解锁高效合同管理系统:优化业务流程,提升管理效率》

随着企业规模的扩大和业务复杂性的增加&#xff0c;合同管理变得愈发重要。合同是企业与客户、供应商、合作伙伴之间的法律约束和商业承诺&#xff0c;而有效的合同管理系统则成为企业提高运营效率、降低风险的关键工具。本文将探讨合同管理系统的重要性以及如何利用合同管理系…

【python】Flask开发感悟

【背景】 做需求做多了&#xff0c;有一个重要的感悟&#xff0c;无关乎技术&#xff0c;但决定着设计思维的方向。 【感悟】 用户往往提出的需求是碎片化的&#xff0c;他们想要的只是一个解决狭窄问题的特定工具。 但是身兼设计师和开发者双重身份的人员不可能看得这么狭窄…

【YashanDB知识库】ycm托管数据库时报错OM host ip:127.0.0.1 is not support join to YCM

问题现象 问题的风险及影响 导致数据库无法托管监控 问题影响的版本 问题发生原因 安装数据库时修改了OM的监听ip为127.0.0.1 解决方法及规避方式 后台修改OM的ip为本机的ip或者0.0.0.0 问题分析和处理过程 1、修改env文件中的om IP地址&#xff0c;修改为0.0.0.0或本机…

阅读完善程序复习(四)2022csp-s

错题&#xff1a; 25. 答案&#xff1a;D 解析&#xff1a;模拟即可。 27. 答案&#xff1a;C 解析&#xff1a;代入一个特例即可&#xff08;基础&#xff0c;可以分析一下排序&#xff09;。 32. 答案&#xff1a;B 解析&#xff1a;模拟即可。 33. 答案&#xff1a;B…

milvus元数据在etcd的存储解析

milvus元数据在etcd的存储解析 数据以key-value形式存在。 大致包含如下一些种类: databasecollectionfieldpartitionindexsegment-indexresource_groupsession database 创建一个数据库会产生2个key&#xff0c;但value是相同的。 key规则: 前缀/root-coord/database/db…

代码随想录算法训练营第六十一天| LeetCode739. 每日温度、496.下一个更大元素 I

一、LeetCode739. 每日温度 题目链接/视频讲解/文章讲解&#xff1a;https://programmercarl.com/0739.%E6%AF%8F%E6%97%A5%E6%B8%A9%E5%BA%A6.html 状态&#xff1a;已解决 1.思路 如果有同学学过单调栈的话这道题还是能看出明显解法的&#xff0c;因为需要寻找任一个元素的右…