当前位置: 首页 > news >正文

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方程(具有纳维尔边界的双调和方程)

我们考虑如下形式的双调和方程的数值解 其中&#xff0c;Ω是欧氏空间中的多边形或多面体域&#xff0c;在其中&#xff0c;d为维度&#xff0c;具有分段利普希茨边界&#xff0c;满足内部锥条件&#xff0c;f(x) ∈ L2(Ω)是给定的函数&#xff0c;∆是标准的拉普拉斯算子。算…...

详细介绍微信小程序app.js

这一节&#xff0c;我们详细介绍app.js 这个文件。这个文件的重要性我就不再赘述&#xff0c;前面已经介绍了。 一、app.js是项目的主控文件 任何一个程序都是需要一个入口的&#xff0c;就好比我们在学c的时候就会有一个main函数&#xff0c;其他语言基本都是一样。很明确的…...

【六 (2)机器学习-EDA探索性数据分析模板】

目录 文章导航一、EDA&#xff1a;二、导入类库三、导入数据四、查看数据类型和缺失情况五、确认目标变量和ID六、查看目标变量分布情况七、特征变量按照数据类型分成定量变量和定性变量八、查看定量变量分布情况九、查看定量变量的离散程度十、查看定量变量与目标变量关系十一…...

Java集合——Map、Set和List总结

文章目录 一、Collection二、Map、Set、List的不同三、List1、ArrayList2、LinkedList 四、Map1、HashMap2、LinkedHashMap3、TreeMap 五、Set 一、Collection Collection 的常用方法 public boolean add(E e)&#xff1a;把给定的对象添加到当前集合中 。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的安装就细不讲了&#xff0c;直接去官网下&#xff0c;实在不会的直接用cdn,省的一番口舌。 cdn.staticfile.net/echarts/4.3.0/echarts.min.js 正入话题哈 什么是EChart&#xff1f; EChart 是一个使用 JavaScript 实现的开源可视化库&#xff0c;Echart支持多种常…...

阿里云8核32G云服务器租用优惠价格表,包括腾讯云和京东云

8核32G云服务器租用优惠价格表&#xff0c;云服务器吧yunfuwuqiba.com整理阿里云8核32G服务器、腾讯云8核32G和京东云8C32G云主机配置报价&#xff0c;腾讯云和京东云是轻量应用服务器&#xff0c;阿里云是云服务器ECS&#xff1a; 阿里云8核32G服务器 阿里云8核32G服务器价格…...

设计模式,工厂方法模式

工厂方法模式概述 工厂方法模式&#xff0c;是对简单工厂模式的进一步抽象和推广。以我个人理解&#xff0c;工厂方法模式就是对生产工厂的抽象&#xff0c;就是用一个生产工厂的工厂来进行目标对象的创建。 工厂方法模式的角色组成和简单工厂方法相比&#xff0c;创建了一个…...

WPF中嵌入3D模型通用结构

背景&#xff1a;wpf本身有提供3D的绘制&#xff0c;但是自己通过代码描绘出3D是比较困难的。3D库helix-toolkit支持调用第三方生成的模型&#xff0c;比如Blender这些&#xff0c;所以在wpf上使用3D就变得非常简单。这里是一个通过helix-toolkit库调用第三方生成的3d模型的样例…...

举个例子说明联邦学习

学习目标&#xff1a; 一周掌握 Java 入门知识 学习内容&#xff1a; 联邦学习是一种机器学习方法&#xff0c;它允许多个参与者协同训练一个共享模型&#xff0c;同时保持各自数据的隐私。 联邦学习概念&#xff08;例子&#xff09;: 假设有三家医院&#xff0c;它们都希望…...

【Python】免费的图片/图标网站

专栏文章索引&#xff1a;Python 有问题可私聊&#xff1a;QQ&#xff1a;3375119339 这里是我收集的几个免费的图片/图标网站&#xff1a; iconfont-阿里巴巴矢量图标库icon&#xff08;.ico&#xff09;INCONFINDER&#xff08;.ico&#xff09;...

Pytorch中的nn.Embedding()

模块的输入是一个索引列表&#xff0c;输出是相应的词嵌入。 Embedding.weight&#xff08;Tensor&#xff09;–形状模块&#xff08;num_embeddings&#xff0c;Embedding_dim&#xff09;的可学习权重&#xff0c;初始化自&#xff08;0&#xff0c;1&#xff09;。 也就是…...

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&#xff08;n<9)&#xff0c;输出由数字组成的直角三角图形。例如&#xff0c;输入5&#xff0c;输出图形如下 nint(input("")) #开始 for i in range(1,n1):for j in range(1,i1):print(j,end"")print()#结束 2.打印字符图形 …...

前端面试题--CSS系列(一)

CSS系列--持续更新中 1.CSS预处理器有哪些类型&#xff0c;有什么区别2.盒模型是什么&#xff0c;有哪两种类型3.css选择器有哪些&#xff0c;优先级是怎样的&#xff0c;哪些属性可以继承4. 说说em/px/rem/vh/vw的区别5.元素实现水平垂直居中的方法有哪些&#xff0c;如果元素…...

VSCode好用插件

由于现在还是使用vue2&#xff0c;所以本文只记录vue2开发中好用的插件。 美化类插件不介绍了&#xff0c;那些貌似对生产力起不到什么大的帮助&#xff0c;纯粹的“唯心主义”罢了&#xff0c;但是如果你有兴趣的话可以查看上一篇博客&#xff1a;VSCode美化 1. vuter 简介&…...

Vue3:对ref、reactive的一个性能优化API

一、情景说明 我们知道&#xff0c;在Vue3中&#xff0c;想要创建响应式的变量&#xff0c;就要用到ref、reactive来包裹一下数据即可。 但是&#xff0c;这里有个损耗性能的地方 就是&#xff0c;被它包裹的数据&#xff0c;都会构建成响应式的&#xff0c;无论多少层次&…...

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 位的地址空间&#xff0c;可以彻底解决 IPv4 地址不足的问题。由于 IPv4 地址是32 位二进制&#xff0c;所能表示的IP 地址个数为 2^32 4 294 967 29640 亿&#xff0c;因而在因特网上约有 40亿个P 地址。 由 32 位的IPv4 升级至 128 位的IPv6&am…...

深入浅出 -- 系统架构之负载均衡Nginx资源压缩

一、Nginx资源压缩 建立在动静分离的基础之上&#xff0c;如果一个静态资源的Size越小&#xff0c;那么自然传输速度会更快&#xff0c;同时也会更节省带宽&#xff0c;因此我们在部署项目时&#xff0c;也可以通过Nginx对于静态资源实现压缩传输&#xff0c;一方面可以节省带宽…...

基于jsp+Spring boot+mybatis的图书管理系统设计和实现

基于jspSpring bootmybatis的图书管理系统设计和实现 博主介绍&#xff1a;多年java开发经验&#xff0c;专注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项目介绍 本项目&#xff08;苍穹外卖&#xff09;是专门为餐饮企业&#xff08;餐厅、饭店&#xff09;定制的一款软件产品&#xff0c;包括 系统管理后台 和 小程序端应用 两部分。其中系统管理后台主要提供给餐饮企业内部员…...

云原生架构(微服务、容器云、DevOps、不可变基础设施、声明式API、Serverless、Service Mesh)

前言 读完本文&#xff0c;你将对云原生下的核心概念微服务、容器云、DevOps、Immutable Infrastructure、Declarative-API、Serverless、Service Mesh 等有一个相对详细的了解&#xff0c;帮助你快速掌握云原生的核心和要点。 因题主资源有限, 这里会选用部分云服务商的组件进…...

基于重写ribbon负载实现灰度发布

项目结构如下 代码如下&#xff1a; pom&#xff1a; <?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*算法 场景题&#xff1a;扫雷如何随机分地雷&#xff0c;怎么安排数字显示 怎么判断一个物体在三角锥内 动作游戏中打击效果怎么处理穿模问题 八叉树了解过吗 骨骼动画…...

idea开发 java web 高校学籍管理系统bootstrap框架web结构java编程计算机网页

一、源码特点 java 高校学籍管理系统是一套完善的完整信息系统&#xff0c;结合java web开发和bootstrap UI框架完成本系统 &#xff0c;对理解JSP java编程开发语言有帮助&#xff0c;系统具有完整的源代码和数据库&#xff0c;系统主要采用B/S模式开发。 前段主要技术 css jq…...

linux之文件系统、inode和动静态库制作和发布

一、背景 1.没有被打开的文件都在磁盘上 --- 磁盘级文件 2.对磁盘级别的文件&#xff0c;我们的侧重点 单个文件角度 -- 这个文件在哪里&#xff0c;有多大&#xff0c;其他属性是什么&#xff1f; 站在系统角度 -- 一共有多少文件&#xff1f;各自属性在哪里&#xff1f…...

C++IO类,输入输出缓冲区,流状态

我们的程序已经使用了很多IO库设施&#xff1a; istream(输入流)类型&#xff0c;提供输入操作。ostream(输出流)类型&#xff0c;提供输出操作。cin&#xff0c;一个istream对象&#xff0c;从标准输入读取数据。写入到标准错误。cout&#xff0c;一个ostream对象&#xff0c…...

机器学习笔记 - 文字转语音技术路线简述以及相关工具不完全清单

一、TTS技术简述 今天的文本到语音转换技术(TTS)的目标已经不仅仅是让机器说话,而是让它们听起来像不同年龄和性别的人类。通常,TTS 系统合成器的质量是从不同方面进行评估的,包括合成语音的清晰度、自然度和偏好,以及人类感知因素,例如可理解性。 1、技术路线 (1)基…...

wordpress自动保存远程图片/seo优化实训总结

MySQL 数据库最佳学习线路脑图&#xff1a; 一、 对MySQL 的认识 认识Mysql数据库 下载安装MySQL软件 在Linux系统环境下安装MySQL MySOL体系结构与存储引擎 MySQL体系结构 Query Cache 详解存储引擎InnoDB体系结构InnoDB的三大特性. 数据库文件 参数文件参数类型错误…...

哪个视频网站做直播销售/泉州关键词排名工具

2019独角兽企业重金招聘Python工程师标准>>> 访问权限修饰符&#xff1a;说明类中成员的访问权限 类中的成员按访问权限分为public&#xff08;公有&#xff09;、protected&#xff08;保护&#xff09;、priavte&#xff08;私有&#xff09;。 公有成员可以被该类…...

深圳哪里有做网站的公司/苏州百度推广代理商

不废话&#xff0c;直接进入主题&#xff0c;本次演示如何启动chrome&#xff0c;且指定默认的profile下载chromedriver地址&#xff1a;http://chromedriver.storage.googleapis.com/index.html?path2.33/特别注意&#xff1a;看notes.txt注意版本号一定要匹配啊&#xff0c;…...

长沙个人做网站/百度识图扫一扫

作者&#xff1a;桂。 链接&#xff1a;http://www.cnblogs.com/xingshansi/articles/6291358.html 声明&#xff1a;转载请注明出处&#xff0c;谢谢。 前言 工欲善其事&#xff0c;必先利其器。在信息爆炸的今天&#xff0c;如何快速有效搜索期望的信息&#xff0c;成了困扰…...

网络服务商基本网络参数/seo基础课程

官方文档的写法&#xff1a; 这里要注意的是&#xff0c;user有两个子路由&#xff0c;分别是profile和posts,如上图设置之后&#xff0c;需要在user.vue里面设置<router-view></router-view>才能显示跳转的子路由页面。...

广告网站设计方案/在哪里可以找到网站

问题描述&#xff1a;当使用Quartus II综合好的文件下载到FPGA时发现&#xff0c;没有识别到下载接口USB_Blaster,如下图所示。 图片中高亮的部分只有显示No Hardware,并没有显示我们需要的USB-Blaster[USB-0]&#xff08;上面截图显示出来USB-Blaster[USB-0]是因为这张图是我…...