当前位置: 首页 > 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;一方面可以节省带宽…...

多云管理“拦路虎”:深入解析网络互联、身份同步与成本可视化的技术复杂度​

一、引言&#xff1a;多云环境的技术复杂性本质​​ 企业采用多云策略已从技术选型升维至生存刚需。当业务系统分散部署在多个云平台时&#xff0c;​​基础设施的技术债呈现指数级积累​​。网络连接、身份认证、成本管理这三大核心挑战相互嵌套&#xff1a;跨云网络构建数据…...

Cursor实现用excel数据填充word模版的方法

cursor主页&#xff1a;https://www.cursor.com/ 任务目标&#xff1a;把excel格式的数据里的单元格&#xff0c;按照某一个固定模版填充到word中 文章目录 注意事项逐步生成程序1. 确定格式2. 调试程序 注意事项 直接给一个excel文件和最终呈现的word文件的示例&#xff0c;…...

调用支付宝接口响应40004 SYSTEM_ERROR问题排查

在对接支付宝API的时候&#xff0c;遇到了一些问题&#xff0c;记录一下排查过程。 Body:{"datadigital_fincloud_generalsaas_face_certify_initialize_response":{"msg":"Business Failed","code":"40004","sub_msg…...

CocosCreator 之 JavaScript/TypeScript和Java的相互交互

引擎版本&#xff1a; 3.8.1 语言&#xff1a; JavaScript/TypeScript、C、Java 环境&#xff1a;Window 参考&#xff1a;Java原生反射机制 您好&#xff0c;我是鹤九日&#xff01; 回顾 在上篇文章中&#xff1a;CocosCreator Android项目接入UnityAds 广告SDK。 我们简单讲…...

sipsak:SIP瑞士军刀!全参数详细教程!Kali Linux教程!

简介 sipsak 是一个面向会话初始协议 (SIP) 应用程序开发人员和管理员的小型命令行工具。它可以用于对 SIP 应用程序和设备进行一些简单的测试。 sipsak 是一款 SIP 压力和诊断实用程序。它通过 sip-uri 向服务器发送 SIP 请求&#xff0c;并检查收到的响应。它以以下模式之一…...

C++课设:简易日历程序(支持传统节假日 + 二十四节气 + 个人纪念日管理)

名人说:路漫漫其修远兮,吾将上下而求索。—— 屈原《离骚》 创作者:Code_流苏(CSDN)(一个喜欢古诗词和编程的Coder😊) 专栏介绍:《编程项目实战》 目录 一、为什么要开发一个日历程序?1. 深入理解时间算法2. 练习面向对象设计3. 学习数据结构应用二、核心算法深度解析…...

[大语言模型]在个人电脑上部署ollama 并进行管理,最后配置AI程序开发助手.

ollama官网: 下载 https://ollama.com/ 安装 查看可以使用的模型 https://ollama.com/search 例如 https://ollama.com/library/deepseek-r1/tags # deepseek-r1:7bollama pull deepseek-r1:7b改token数量为409622 16384 ollama命令说明 ollama serve #&#xff1a…...

R 语言科研绘图第 55 期 --- 网络图-聚类

在发表科研论文的过程中&#xff0c;科研绘图是必不可少的&#xff0c;一张好看的图形会是文章很大的加分项。 为了便于使用&#xff0c;本系列文章介绍的所有绘图都已收录到了 sciRplot 项目中&#xff0c;获取方式&#xff1a; R 语言科研绘图模板 --- sciRplothttps://mp.…...

C语言中提供的第三方库之哈希表实现

一. 简介 前面一篇文章简单学习了C语言中第三方库&#xff08;uthash库&#xff09;提供对哈希表的操作&#xff0c;文章如下&#xff1a; C语言中提供的第三方库uthash常用接口-CSDN博客 本文简单学习一下第三方库 uthash库对哈希表的操作。 二. uthash库哈希表操作示例 u…...

【C++】纯虚函数类外可以写实现吗?

1. 答案 先说答案&#xff0c;可以。 2.代码测试 .h头文件 #include <iostream> #include <string>// 抽象基类 class AbstractBase { public:AbstractBase() default;virtual ~AbstractBase() default; // 默认析构函数public:virtual int PureVirtualFunct…...