matlab:五点中心差分求解Navier边界的Biharmonic方程(具有纳维尔边界的双调和方程)
我们考虑如下形式的双调和方程的数值解
其中,Ω是欧氏空间中的多边形或多面体域,在其中,d为维度,具有分段利普希茨边界,满足内部锥条件,f(x) ∈ L2(Ω)是给定的函数,∆是标准的拉普拉斯算子。算子∆u(x)和∆2u(x)表示为
巧妙地将双调和方程(1.1)分解为两个Possion方程,传统的数值方法如有限元法(FEM)和有限差分法(FDM)在计算资源和时间复杂度较小的情况下表现良好。通过引入辅助变量w = −∆u,可以将四阶方程(1.1)重写为一对二阶方程:
或者引入变量w = ∆u,得到
那么,我们的目标为寻找一对函数(w,u),而不是找到原始问题(1.1)的解。如下我们以g=0和h=0为例,利用五点中心差分求解上面的双调和方程。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Matrix method for Biharmonic Equation %%%%
%%% u_{xxxx} + u_{yyyy} + 2*u_{xx}*u_{yy} = f(x, y) %%%%
%%% Omega = 0 < x < 1, 0 < y < 1 %%%%
%%% u(x, y) = 0 on boundary, %%%%
%%% Exact soln: u(x, y) = sin(pi*x)*sin(pi*y) %%%%
%%% Here f(x, y) = 4*pi^4*sin(pi*x)*sin(pi*y); %%%%
%%% Course: Xi'an Li on 12.08 2023 %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear all
clc
close allftsz = 20;x_l = -1.0;
x_r = 1.0;
y_b = -1.0;
y_t = 1.0;q = 6;
Num = 2^q+1;
NNN = Num*Num;point_num2x = Num;
point_num2y = Num;fside = @(x, y) 4*pi^4*sin(pi*x).*sin(pi*y);
utrue = @(x, y) sin(pi*x).*sin(pi*y);hx = (x_r-x_l)/point_num2x;
X = zeros(point_num2x-1,1);
for i=1:point_num2x-1X(i) = x_l+i*hx;
endhy = (y_t-y_b)/point_num2y;
Y=zeros(point_num2y-1,1);
for i=1:point_num2y-1Y(i) = y_b+i*hy;
end
[meshX, meshY] = meshgrid(X, Y);tic;
Unumberi = FDM2Biharmonic_Couple2Navier_Zero(point_num2x, point_num2y,...x_l, x_r, y_b, y_t, fside);
fprintf('%s%s%s\n','运行时间:',toc,'s')
U_exact = utrue(meshX, meshY);
meshErr = abs(U_exact - Unumberi);rel_err = sum(sum(meshErr))/sum(sum(abs(U_exact)));
fprintf('%s%s\n','相对误差:',rel_err)figure('name','Exact')
axis tight;
h = surf(meshX, meshY, U_exact','Edgecolor', 'none');
hold on
title('Exact Solu')
% xlabel('$x$', 'Fontsize', 20, 'Interpreter', 'latex')
% ylabel('$y$', 'Fontsize', 20, 'Interpreter', 'latex')
% zlabel('$Error$', 'Fontsize', 20, 'Interpreter', 'latex')
hold on
set(gca, 'XMinortick', 'off', 'YMinorTick', 'off', 'Fontsize', ftsz);
set(gcf, 'Renderer', 'zbuffer');
hold on
% colorbar;
% caxis([0 0.00012])
hold onfigure('name','Absolute Error')
axis tight;
h = surf(meshX, meshY, meshErr','Edgecolor', 'none');
hold on
title('Absolute Error')
% xlabel('$x$', 'Fontsize', 20, 'Interpreter', 'latex')
% ylabel('$y$', 'Fontsize', 20, 'Interpreter', 'latex')
% zlabel('$Error$', 'Fontsize', 20, 'Interpreter', 'latex')
hold on
set(gca, 'XMinortick', 'off', 'YMinorTick', 'off', 'Fontsize', ftsz);
set(gcf, 'Renderer', 'zbuffer');
hold on
% colorbar;
% caxis([0 0.00012])
hold onif q==6txt2result = 'result2fdm_mesh6.txt';
elseif q==7txt2result = 'result2fdm_mesh7.txt';
elseif q==8txt2result = 'result2fdm_mesh8.txt';
elseif q==9txt2result = 'result2fdm_mesh9.txt';
endfop = fopen(txt2result, 'wt');fprintf(fop,'%s%s%s\n','运行时间:',toc,'s');
fprintf(fop,'%s%d\n','内部网格点数目:',Num-1);
fprintf(fop,'%s%s\n','相对误差:',rel_err);
被调用的求解函数如下:
function Uapp = FDM2Biharmonic_Couple2Navier_Zero(Nx, Ny, xleft, xright, ybottom, ytop, fside)format long;% Define the step sizes and create the grid without boundary pointshx = (xright-xleft)/Nx; x = zeros(Nx-1,1);for ix=1:Nx-1x(ix) = xleft+ix*hx;endhy = (ytop-ybottom)/Ny; y=zeros(Ny-1,1);for iy=1:Ny-1y(iy) = ybottom+iy*hy;end% Define the source termsource_term = fside;% Initialize the coefficient matrix A and the right-hand side vector FN = (Ny-1)*(Nx-1);A = sparse(N,N); FV = zeros(N,1);% Loop through each inner grid point, Apply finite difference scheme (central differences)hx1 = hx*hx; hy1 = hy*hy; for jv = 1:Ny-1for iv=1:Nx-1kv = iv + (jv-1)*(Nx-1);FV(kv) = fside(x(iv),y(jv));A(kv,kv) = -2/hx1 -2/hy1;%-- x direction --------------if iv == 1A(kv,kv+1) = 1/hx1;elseif iv==Nx-1A(kv,kv-1) = 1/hx1;elseA(kv,kv-1) = 1/hx1;A(kv,kv+1) = 1/hx1;endend%-- y direction --------------if jv == 1A(kv,kv+Nx-1) = 1/hy1;elseif jv==Ny-1A(kv,kv-(Nx-1)) = 1/hy1;elseA(kv,kv-(Nx-1)) = 1/hy1;A(kv,kv+Nx-1) = 1/hy1;endendendendV = A\FV;B = sparse(N,N); FU = zeros(N,1);% Loop through each inner grid point, Apply finite difference scheme (central differences)for ju = 1:Ny-1for iu=1:Nx-1ku = iu + (ju-1)*(Nx-1);FV(ku) = V(ku);B(ku,ku) = -2/hx1 -2/hy1;%-- x direction --------------if iu == 1B(ku,ku+1) = 1/hx1;elseif iu==Nx-1B(ku,ku-1) = 1/hx1;elseB(ku,ku-1) = 1/hx1;B(ku,ku+1) = 1/hx1;endend%-- y direction --------------if ju == 1B(ku,ku+Nx-1) = 1/hy1;elseif ju==Ny-1B(ku,ku-(Nx-1)) = 1/hy1;elseB(ku,ku-(Nx-1)) = 1/hy1;B(ku,ku+Nx-1) = 1/hy1;endendendendU = B\FV;%--- Transform back to (i,j) form to plot the solution ---j = 1;for k=1:Ni = k - (j-1)*(Nx-1) ;Uapp(i,j) = U(k);j = fix(k/(Nx-1)) + 1;end
end
结果如下:
运行时间:6.574370e-02s
相对误差:1.558669e-03
相关文章:
matlab:五点中心差分求解Navier边界的Biharmonic方程(具有纳维尔边界的双调和方程)
我们考虑如下形式的双调和方程的数值解 其中,Ω是欧氏空间中的多边形或多面体域,在其中,d为维度,具有分段利普希茨边界,满足内部锥条件,f(x) ∈ L2(Ω)是给定的函数,∆是标准的拉普拉斯算子。算…...
详细介绍微信小程序app.js
这一节,我们详细介绍app.js 这个文件。这个文件的重要性我就不再赘述,前面已经介绍了。 一、app.js是项目的主控文件 任何一个程序都是需要一个入口的,就好比我们在学c的时候就会有一个main函数,其他语言基本都是一样。很明确的…...
【六 (2)机器学习-EDA探索性数据分析模板】
目录 文章导航一、EDA:二、导入类库三、导入数据四、查看数据类型和缺失情况五、确认目标变量和ID六、查看目标变量分布情况七、特征变量按照数据类型分成定量变量和定性变量八、查看定量变量分布情况九、查看定量变量的离散程度十、查看定量变量与目标变量关系十一…...
Java集合——Map、Set和List总结
文章目录 一、Collection二、Map、Set、List的不同三、List1、ArrayList2、LinkedList 四、Map1、HashMap2、LinkedHashMap3、TreeMap 五、Set 一、Collection Collection 的常用方法 public boolean add(E e):把给定的对象添加到当前集合中 。public void clear(…...
Python TensorFlow 2.6 获取 MNIST 数据
Python TensorFlow 2.6 获取 MNIST 数据 2 Python TensorFlow 2.6 获取 MNIST 数据1.1 获取 MNIST 数据1.2 检查 MNIST 数据 2 Python 将npz数据保存为txt3 Java 获取数据并使用SVM训练4 Python 测试SVM准确度 2 Python TensorFlow 2.6 获取 MNIST 数据 1.1 获取 MNIST 数据 …...
EChart简单入门
echart的安装就细不讲了,直接去官网下,实在不会的直接用cdn,省的一番口舌。 cdn.staticfile.net/echarts/4.3.0/echarts.min.js 正入话题哈 什么是EChart? EChart 是一个使用 JavaScript 实现的开源可视化库,Echart支持多种常…...
阿里云8核32G云服务器租用优惠价格表,包括腾讯云和京东云
8核32G云服务器租用优惠价格表,云服务器吧yunfuwuqiba.com整理阿里云8核32G服务器、腾讯云8核32G和京东云8C32G云主机配置报价,腾讯云和京东云是轻量应用服务器,阿里云是云服务器ECS: 阿里云8核32G服务器 阿里云8核32G服务器价格…...
设计模式,工厂方法模式
工厂方法模式概述 工厂方法模式,是对简单工厂模式的进一步抽象和推广。以我个人理解,工厂方法模式就是对生产工厂的抽象,就是用一个生产工厂的工厂来进行目标对象的创建。 工厂方法模式的角色组成和简单工厂方法相比,创建了一个…...
WPF中嵌入3D模型通用结构
背景:wpf本身有提供3D的绘制,但是自己通过代码描绘出3D是比较困难的。3D库helix-toolkit支持调用第三方生成的模型,比如Blender这些,所以在wpf上使用3D就变得非常简单。这里是一个通过helix-toolkit库调用第三方生成的3d模型的样例…...
举个例子说明联邦学习
学习目标: 一周掌握 Java 入门知识 学习内容: 联邦学习是一种机器学习方法,它允许多个参与者协同训练一个共享模型,同时保持各自数据的隐私。 联邦学习概念(例子): 假设有三家医院,它们都希望…...
【Python】免费的图片/图标网站
专栏文章索引:Python 有问题可私聊:QQ:3375119339 这里是我收集的几个免费的图片/图标网站: iconfont-阿里巴巴矢量图标库icon(.ico)INCONFINDER(.ico)...
Pytorch中的nn.Embedding()
模块的输入是一个索引列表,输出是相应的词嵌入。 Embedding.weight(Tensor)–形状模块(num_embeddings,Embedding_dim)的可学习权重,初始化自(0,1)。 也就是…...
WebSocketServer后端配置,精简版
首先需要maven配置 <dependency><groupId>org.springframework.boot</groupId><artifactId>spring-boot-starter-websocket</artifactId><version>2.1.3.RELEASE</version></dependency> 然后加上配置类 这段代码是一个Spri…...
Python程序设计 多重循环(二)
1.打印数字图形 输入n(n<9),输出由数字组成的直角三角图形。例如,输入5,输出图形如下 nint(input("")) #开始 for i in range(1,n1):for j in range(1,i1):print(j,end"")print()#结束 2.打印字符图形 …...
前端面试题--CSS系列(一)
CSS系列--持续更新中 1.CSS预处理器有哪些类型,有什么区别2.盒模型是什么,有哪两种类型3.css选择器有哪些,优先级是怎样的,哪些属性可以继承4. 说说em/px/rem/vh/vw的区别5.元素实现水平垂直居中的方法有哪些,如果元素…...
VSCode好用插件
由于现在还是使用vue2,所以本文只记录vue2开发中好用的插件。 美化类插件不介绍了,那些貌似对生产力起不到什么大的帮助,纯粹的“唯心主义”罢了,但是如果你有兴趣的话可以查看上一篇博客:VSCode美化 1. vuter 简介&…...
Vue3:对ref、reactive的一个性能优化API
一、情景说明 我们知道,在Vue3中,想要创建响应式的变量,就要用到ref、reactive来包裹一下数据即可。 但是,这里有个损耗性能的地方 就是,被它包裹的数据,都会构建成响应式的,无论多少层次&…...
Python 用pygame简简单单实现一个打砖块
# -*- coding: utf-8 -*- # # # Copyright (C) 2024 , Inc. All Rights Reserved # # # Time : 2024/3/30 14:34 # Author : 赫凯 # Email : hekaiiii163.com # File : ballgame.py # Software: PyCharm import math import randomimport pygame import sys#…...
软考113-上午题-【计算机网络】-IPv6、无线网络、Windows命令
一、IPv6 IPv6 具有长达 128 位的地址空间,可以彻底解决 IPv4 地址不足的问题。由于 IPv4 地址是32 位二进制,所能表示的IP 地址个数为 2^32 4 294 967 29640 亿,因而在因特网上约有 40亿个P 地址。 由 32 位的IPv4 升级至 128 位的IPv6&am…...
深入浅出 -- 系统架构之负载均衡Nginx资源压缩
一、Nginx资源压缩 建立在动静分离的基础之上,如果一个静态资源的Size越小,那么自然传输速度会更快,同时也会更节省带宽,因此我们在部署项目时,也可以通过Nginx对于静态资源实现压缩传输,一方面可以节省带宽…...
基于jsp+Spring boot+mybatis的图书管理系统设计和实现
基于jspSpring bootmybatis的图书管理系统设计和实现 博主介绍:多年java开发经验,专注Java开发、定制、远程、文档编写指导等,csdn特邀作者、专注于Java技术领域 作者主页 央顺技术团队 Java毕设项目精品实战案例《1000套》 欢迎点赞 收藏 ⭐留言 文末获…...
Pytorch转onnx
pytorch 转 onnx 模型需要函数 torch.onnx.export。 def export(model: Union[torch.nn.Module, torch.jit.ScriptModule, torch.jit.ScriptFunction],args: Union[Tuple[Any, ...], torch.Tensor],f: Union[str, io.BytesIO],export_params: bool True,verbose: bool False…...
苍穹外卖——项目搭建
一、项目介绍以及环境搭建 1.苍穹外卖项目介绍 1.1项目介绍 本项目(苍穹外卖)是专门为餐饮企业(餐厅、饭店)定制的一款软件产品,包括 系统管理后台 和 小程序端应用 两部分。其中系统管理后台主要提供给餐饮企业内部员…...
云原生架构(微服务、容器云、DevOps、不可变基础设施、声明式API、Serverless、Service Mesh)
前言 读完本文,你将对云原生下的核心概念微服务、容器云、DevOps、Immutable Infrastructure、Declarative-API、Serverless、Service Mesh 等有一个相对详细的了解,帮助你快速掌握云原生的核心和要点。 因题主资源有限, 这里会选用部分云服务商的组件进…...
基于重写ribbon负载实现灰度发布
项目结构如下 代码如下: pom: <?xml version"1.0" encoding"UTF-8"?> <project xmlns"http://maven.apache.org/POM/4.0.0"xmlns:xsi"http://www.w3.org/2001/XMLSchema-instance"xsi:schemaLocat…...
无端科技一面(生死狙击项目组 战斗客户端 40min)
自我介绍 实习经历询问 项目询问 TCP和UDP的区别 什么情况会用到UDP 大小端 寻路算法了解多少 A*算法 场景题:扫雷如何随机分地雷,怎么安排数字显示 怎么判断一个物体在三角锥内 动作游戏中打击效果怎么处理穿模问题 八叉树了解过吗 骨骼动画…...
idea开发 java web 高校学籍管理系统bootstrap框架web结构java编程计算机网页
一、源码特点 java 高校学籍管理系统是一套完善的完整信息系统,结合java web开发和bootstrap UI框架完成本系统 ,对理解JSP java编程开发语言有帮助,系统具有完整的源代码和数据库,系统主要采用B/S模式开发。 前段主要技术 css jq…...
linux之文件系统、inode和动静态库制作和发布
一、背景 1.没有被打开的文件都在磁盘上 --- 磁盘级文件 2.对磁盘级别的文件,我们的侧重点 单个文件角度 -- 这个文件在哪里,有多大,其他属性是什么? 站在系统角度 -- 一共有多少文件?各自属性在哪里?…...
C++IO类,输入输出缓冲区,流状态
我们的程序已经使用了很多IO库设施: istream(输入流)类型,提供输入操作。ostream(输出流)类型,提供输出操作。cin,一个istream对象,从标准输入读取数据。写入到标准错误。cout,一个ostream对象,…...
机器学习笔记 - 文字转语音技术路线简述以及相关工具不完全清单
一、TTS技术简述 今天的文本到语音转换技术(TTS)的目标已经不仅仅是让机器说话,而是让它们听起来像不同年龄和性别的人类。通常,TTS 系统合成器的质量是从不同方面进行评估的,包括合成语音的清晰度、自然度和偏好,以及人类感知因素,例如可理解性。 1、技术路线 (1)基…...
wordpress自动保存远程图片/seo优化实训总结
MySQL 数据库最佳学习线路脑图: 一、 对MySQL 的认识 认识Mysql数据库 下载安装MySQL软件 在Linux系统环境下安装MySQL MySOL体系结构与存储引擎 MySQL体系结构 Query Cache 详解存储引擎InnoDB体系结构InnoDB的三大特性. 数据库文件 参数文件参数类型错误…...
哪个视频网站做直播销售/泉州关键词排名工具
2019独角兽企业重金招聘Python工程师标准>>> 访问权限修饰符:说明类中成员的访问权限 类中的成员按访问权限分为public(公有)、protected(保护)、priavte(私有)。 公有成员可以被该类…...
深圳哪里有做网站的公司/苏州百度推广代理商
不废话,直接进入主题,本次演示如何启动chrome,且指定默认的profile下载chromedriver地址:http://chromedriver.storage.googleapis.com/index.html?path2.33/特别注意:看notes.txt注意版本号一定要匹配啊,…...
长沙个人做网站/百度识图扫一扫
作者:桂。 链接:http://www.cnblogs.com/xingshansi/articles/6291358.html 声明:转载请注明出处,谢谢。 前言 工欲善其事,必先利其器。在信息爆炸的今天,如何快速有效搜索期望的信息,成了困扰…...
网络服务商基本网络参数/seo基础课程
官方文档的写法: 这里要注意的是,user有两个子路由,分别是profile和posts,如上图设置之后,需要在user.vue里面设置<router-view></router-view>才能显示跳转的子路由页面。...
广告网站设计方案/在哪里可以找到网站
问题描述:当使用Quartus II综合好的文件下载到FPGA时发现,没有识别到下载接口USB_Blaster,如下图所示。 图片中高亮的部分只有显示No Hardware,并没有显示我们需要的USB-Blaster[USB-0](上面截图显示出来USB-Blaster[USB-0]是因为这张图是我…...