Cranck-Nicolson隐式方法解线性双曲型方程

news2024/12/23 22:54:59

Cranck-Nicolson隐式方法解线性双曲型方程
Cranck-Nicolson方法在抛物型方程里面比较常用,双曲型方程例子不多,该方法是二阶精度,无条件稳定,然而,数值震荡比较明显,特别是时间演化比较大以及courant数比较大的时候。
对于这类方程的隐式解法,系数矩阵是三对角矩阵,每次固定时间t,通过解方程的方法解出来下一个时间步上的函数值,比如X方向分了N个格点,去掉边界条件,每个时间步N-2个未知数,构成三对角矩阵。如果采用追赶法求解,courant数需要小于2
网上的例子不多,我这里写了一下,供大家参考
在这里插入图片描述
CN方法离散成差分方程组,注意对于边界条件旁边的未知数边界项需要移动到等式右边
在这里插入图片描述

在这里插入图片描述在这里插入图片描述
这里外边界直接采用精确解作为边界条件

import matplotlib
import math
matplotlib.use('TkAgg')
import numpy as np  
import matplotlib.pyplot as plt


def catchup(n,a0,b0,c0,F0):             #输入三个对角线上的值a,b,c,F0是等号右边的向量
#注意,追赶法要求:|b1|>|c1|,|bn|>|an|, |bi|>|ai|+|ci|

    m = np.zeros(F0.size, dtype=float)#中间变量
    xx = np.zeros(F0.size, dtype=float)
    for i in range(1, n): #正序,追
        m[i] = a0[i]/b0[i-1]#          a0是左副对角线,c0是右副对角线

        b0[i] = np.copy(b0[i]) - m[i]*c0[i-1]#b是主对角线
        F0[i] = np.copy(F0[i]) - m[i]*F0[i-1]
    xx[-1] = F0[-1]/b0[-1]
    for j in range (n-2, -1, -1):#倒序,赶
        xx[j] = (F0[j]-c0[j]*xx[j+1])/b0[j]
    return xx


def CN(x,t,a0,u):#Cranck-Nicolson格式
    n = (x.size-2)#  x方向未知数数量,x方向有两个边界作为已知条件,因此未知数是格点数-2
    m = t.size-1   #t方向未知数量,t方向未知数-1

    courant = a0*(t[2]-t[1])/(x[2]-x[1])#courant数,虽然CN算法对courant数没有限制,但追赶法有限制
    if courant > 2:
        print('courant number is too large',courant )
    #由于追赶法要求三条对角线|bi|>|ai|+|ci|,对应从courant<2
    print('courant=',courant)

    b = np.zeros(n, dtype=float)#方程组的系数矩阵,b是主对角线,ac是两条副对角线,这是一个三对角矩阵
    a = np.zeros(n, dtype=float)
    c = np.zeros(n, dtype=float)
    f = np.zeros(n, dtype=float)#差分方程等号的右边,和对角线上元素数量一致

    for j in range (0,m):

        for i in range(0,n): #对三条对角线分别赋值,对方程组等号右边f也赋值,注意第一行和最后一行需要额外+-courant/4.*u[m+1,0]
            b[i] = 1.#     主对角线
            if i >0:   
                a[i] = -1.*courant/4.
            if i<(n-1):
                c[i] =courant/4.
            f[i] = u[j,i+1]-courant/4.*(u[j,i+2]-u[j,i])

        '''#这里是系数矩阵,实际上不需要,可以打印出来看一下
        A = np.zeros((n,n),dtype=float)
        for i in range(0,n):         
            A[i,i] = 1.#     主对角线
            if i >0:
                A[i-1,i] = -1.*courant/2.
            if i<(n-1):
                A[i,i+1] =courant/2.
            f[i] = u[j,i+1]-courant/4.*(u[j,i+2]-u[j,i])
            f[0] = np.copy(f[0])+courant/4.*u[j+1,0]
            f[-1] = np.copy(f[-1])-courant/4.*u[j+1,-1]
        aa,bb,cc,=get_tridiagonal2(A,f)
        dd=np.copy(f)
        '''

        #u[m+1,-1] = u[m,-1]#+((t[2]-t[1])/(x[2]-x[1]))*(u[m,-1]-u[m,-2])
        f[0] =  u[j,1]-courant/4.*(u[j,2]-u[j,0])+courant/4.*u[j+1,0]#对方程组等号右边f首元素单独赋值(见公式)
        f[-1] = np.copy(f[-1])-courant/4.*u[j+1,-1]#对方程组等号右边f末元素单独赋值(见公式)
        
        result = catchup(n, a, b, c, f)#采用追赶法,只需输入三条对角线和方程右边
        u[j+1, 1:(result.size+1)] = np.copy(result)
    return u

def plot(x, t, result):
    plt.figure()
    plt.plot(x, result[-1,:])
    plt.show()
    plt.close()

    plt.figure()
    plt.contourf(x, t, result, 50, cmap = 'jet')
    plt.colorbar()
    plt.savefig('CN.png')
    plt.show()
    plt.close()
    return 0


#初始化
a0 = 250.
x = np.linspace(0., 400., 200)
t = np.linspace(0., 0.5, 200)

u = np.zeros((t.size,x.size), dtype=float)
u[0,:] = 100.*np.cos(math.pi*x/60.)
u[:,0] = 100.*np.cos(math.pi*(0.- a0*t)/60.)#100.*np.cos(-1.*math.pi*a0*t/60.)
u[:,-1] = 100.*np.cos(math.pi*(x[-1] - a0*t)/60.)

result = CN(x, t, a0, u)

plot(x, t, result)

在这里插入图片描述

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

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

相关文章

基于Spring Boot的实验室管理系统设计与实现

基于Spring Boot的实验室管理系统设计与实现 开发语言&#xff1a;Java框架&#xff1a;springbootJDK版本&#xff1a;JDK1.8数据库工具&#xff1a;Navicat11开发软件&#xff1a;eclipse/myeclipse/idea 系统部分展示 管理员登录界面&#xff0c;管理员通过后台登录窗口进…

Maven介绍 主要包括Maven的基本介绍,作用,以及对应的Maven模型,可以对Maven有一个基本的了解

1、Maven介绍 1.1 什么是Maven Maven是Apache旗下的一个开源项目&#xff0c;是一款用于管理和构建java项目的工具。 官网&#xff1a;https://maven.apache.org/ Apache 软件基金会&#xff0c;成立于1999年7月&#xff0c;是目前世界上最大的最受欢迎的开源软件基金会&…

Visual Studio Code基础:打开一个编辑器(文件)时,覆盖了原编辑器

相关阅读 VS codehttps://blog.csdn.net/weixin_45791458/category_12658212.html?spm1001.2014.3001.5482 在使用vscode时&#xff0c;偶尔会出现这样的问题&#xff1a;打开了某个编辑器&#xff08;文件&#xff0c;下面统称文件&#xff09;后&#xff0c;再打开其他文件…

【触摸案例-控件不能响应的情况 Objective-C语言】

一、接下来,我们来说这个“控件不能响应的情况”, 1.素材里边,有一个“不接受用户交互的情况”,这么一个代码,把它打开, 把这个项目啊,复制过来,改一个名字,叫做“04-控件不能响应的情况”, 打开之后,command + R,运行一下, 在storyboard上,你也可以看得出来,我…

系统服务(22年国赛)—— nmcli命令部署VXLAN

前言&#xff1a;原文在我的博客网站中&#xff0c;持续更新数通、系统方面的知识&#xff0c;欢迎来访&#xff01; 系统服务&#xff08;22年国赛&#xff09;—— VXLAN服务部署https://myweb.myskillstree.cn/118.html 目录 题目&#xff1a; AppSrv 关闭防火墙和SEli…

哈夫曼编码---一种无损数据压缩算法

哈夫曼编码是一种无损数据压缩算法&#xff0c;该算法在数据压缩&#xff0c;存储和网络传输等领域广泛引用&#xff0c;对互联网的发展也产生了深远的影响。 大家熟知的数据无损压缩软件&#xff0c;如WinRAR&#xff0c;gzip&#xff0c;bzip&#xff0c;lzw&#xff0c;7-z…

ThreeJs 环境配置及遇到问题的解决方法

一、环境搭建 ThreeJs在实际在实际使用中更多的是结合框架开发例如&#xff1a;vue框架、react框架&#xff0c;在使用时需要配置开发环境&#xff0c;本文使用的是vscode ThreeJs NodeJs vue 1、ThreeJs安装 下载路径&#xff1a;GitHub - mrdoob/three.js: JavaScript…

如何进行制造设备数据汇集,发挥数据的价值?

数字化转型正深刻推动制造企业实现远程监控、提高生产效率、降低生产成本、优化产品质量及明晰精细化方向。并且工业互联网的发展离不开工业数据的应用&#xff0c;而制造设备数据汇集正是应用的基础。但制造设备数据汇集存在以下难点及痛点&#xff1a; 1、安全把控难 关键的…

windows日志怎么打开/查看?

windows日志里可以查看到系统进行的各种操作,包括正常开关机记录、dhcp配置警告信息等等,不过很多小伙伴并不知道怎么打开windows日志,为此为大家整理了三种快速打开windows系统日志的方法,大家有需要的话赶紧来看看吧。 目录 一、方法1 二、方法2 三、方法3 (推荐) 一…

基于SpringBoot+Vue乡村养老服务管理系统

项目介绍&#xff1a; 使用旧方法对乡村养老服务管理系统登录的信息进行系统化管理已经不再让人们信赖了&#xff0c;把现在的网络信息技术运用在乡村养老服务管理系统登录的管理上面可以解决许多信息管理上面的难题&#xff0c;比如处理数据时间很长&#xff0c;数据存在错误…

如何用Python语言实现远程控制4路控制器/断路器

如何用Python语言实现远程控制4路控制器/断路器呢&#xff1f; 本文描述了使用Python语言调用HTTP接口&#xff0c;实现控制4路控制器/断路器&#xff0c;支持4路输出&#xff0c;均可独立控制&#xff0c;可接入各种电器。 可选用产品&#xff1a;可根据实际场景需求&#xf…

网红大佬的面子,高阶智驾的里子 | 2024北京车展

相关阅读&#xff1a;2023北京车展 《没有争奇斗艳的车模&#xff0c;只有往死里卷的智能汽车》。 文&#xff5c;刘俊宏 李想、李斌绑定“车圈新顶流”雷军互相抬轿子&#xff0c;红衣大叔周鸿祎高情商点评各家汽车新品...... 为了流量&#xff0c;今年车企大佬们比任何时候…

2024李卜常识王小晨申论类比刷题课

2024年&#xff0c;李卜常识与王小晨申论类比刷题课成为备考公务员考试的热门选择。李卜老师以其深厚的学识&#xff0c;为学员们剖析常识的精髓&#xff1b;而王小晨老师则通过类比刷题的方式&#xff0c;帮助学员们掌握申论的技巧。这两门课程相互补充&#xff0c;让学员们在…

「51媒体」政企宣传邀请媒体的作用?

传媒如春雨&#xff0c;润物细无声&#xff0c;大家好&#xff0c;我是51媒体网胡老师。 政企宣传邀请媒体的作用主要体现在以下几个方面&#xff1a; 提升品牌知名度&#xff1a;通过媒体广泛报道活动内容、亮点及企业形象&#xff0c;可以提升企业或政府的品牌知名度。 增加…

kaggle(4) Regression with an Abalone Dataset 鲍鱼数据集的回归

kaggle&#xff08;4&#xff09; Regression with an Abalone Dataset 鲍鱼数据集的回归 import pandas as pd import numpy as npimport xgboost import lightgbm import optuna import catboostfrom sklearn.model_selection import train_test_split from sklearn.metrics …

JS----前端不同格式的 UUID生成

前端 UUID生成 export const generateUUID () > {let timestamp new Date().getTime()let performanceNow (typeof performance ! undefined && performance.now && performance.now() * 1000) || 0return xxxxxxxx-xxxx-4xxx-yxxx-xxxxxxxxxxxx.replac…

【Linux】进程的控制①之进程创建与进程退出

一 、进程的创建 1、fork函数 fork函数功能&#xff1a;从已经存在的进程中创建一个新进程。新进程为子进程&#xff0c;原进程为父进程。 fork函数创建进程过后&#xff0c;父子进程代码和数据是共享的。在前面也讲过。 2.函数的返回值 如果进程创建成功&#xff0c;给父进…

ResNeXt网络结构

一、简介 在ResNet的基础上&#xff0c;对残差结构的block进行了更新。 ResNeXt网络是一种深度神经网络架构&#xff0c;可以视为对ResNet&#xff08;残差网络&#xff09;的改进和升级。ResNeXt结合了VGG网络的堆叠相同基础模块的策略以及Inception系列网络中的split-trans…

软件模型(简洁明了)

《 软件测试基础持续更新中》 一、软件开发模型 1.1 大爆炸模型 优点&#xff1a;思路简单&#xff0c; 通常可能是开发者的“突发奇 想” 缺点&#xff1a;开发过程是非工程化的&#xff0c;随意性大&#xff0c;结果不可预知 测试&#xff1a;开发任务完成后&#xff0c;…

Docker学习(二十五)构建 Arthas 基础镜像

目录 一、简介二、构建基础镜像2.1 下载 Arthas2.2 编写 Dockerfile2.3 构建镜像2.4 创建容器2.5 测试 一、简介 Arthas 是一款由 阿里巴巴 开发的 线上监控诊断工具。通过全局视角实时查看应用负载、内存、GC、线程等信息&#xff0c;能在不修改代码的情况下&#xff0c;对业…