matlab:五点中心差分求解Navier边界的Biharmonic方程(具有纳维尔边界的双调和方程)

本文主要是介绍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方程(具有纳维尔边界的双调和方程)的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/881880

相关文章

跨国公司撤出在华研发中心的启示:中国IT产业的挑战与机遇

近日,IBM中国宣布撤出在华的两大研发中心,这一决定在IT行业引发了广泛的讨论和关注。跨国公司在华研发中心的撤出,不仅对众多IT从业者的职业发展带来了直接的冲击,也引发了人们对全球化背景下中国IT产业竞争力和未来发展方向的深思。面对这一突如其来的变化,我们应如何看待跨国公司的决策?中国IT人才又该如何应对?中国IT产业将何去何从?本文将围绕这些问题展开探讨。 跨国公司撤出的背景与

poj 3159 (spfa差分约束最短路) poj 1201

poj 3159: 题意: 每次给出b比a多不多于c个糖果,求n最多比1多多少个糖果。 解析: 差分约束。 这个博客讲差分约束讲的比较好: http://www.cnblogs.com/void/archive/2011/08/26/2153928.html 套个spfa。 代码: #include <iostream>#include <cstdio>#i

【生成模型系列(初级)】嵌入(Embedding)方程——自然语言处理的数学灵魂【通俗理解】

【通俗理解】嵌入(Embedding)方程——自然语言处理的数学灵魂 关键词提炼 #嵌入方程 #自然语言处理 #词向量 #机器学习 #神经网络 #向量空间模型 #Siri #Google翻译 #AlexNet 第一节:嵌入方程的类比与核心概念【尽可能通俗】 嵌入方程可以被看作是自然语言处理中的“翻译机”,它将文本中的单词或短语转换成计算机能够理解的数学形式,即向量。 正如翻译机将一种语言

poj 3169 spfa 差分约束

题意: 给n只牛,这些牛有些关系。 ml个关系:fr 与 to 牛间的距离要小于等于 cost。 md个关系:fr 与 to 牛间的距离要大于等于 cost。 隐含关系: d[ i ] <= d[ i + 1 ] 解析: 用以上关系建图,求1-n间最短路即可。 新学了一种建图的方法。。。。。。 代码: #include <iostream>#include

POJ 1364差分约束

给出n个变量,m个约束公式 Sa + Sa+1 + .... + Sa+b < ki or > ki ,叫你判断是否存在着解满足这m组约束公式。 Sa + Sa+1   +   .+ Sa+b =  Sum[a+b] - Sum[a-1]  . 注意加入源点n+1 。 public class Main {public static void main(Strin

matlab读取NC文件(含group)

matlab读取NC文件(含group): NC文件数据结构: 代码: % 打开 NetCDF 文件filename = 'your_file.nc'; % 替换为你的文件名% 使用 netcdf.open 函数打开文件ncid = netcdf.open(filename, 'NC_NOWRITE');% 查看文件中的组% 假设我们想读取名为 "group1" 的组groupName

利用matlab bar函数绘制较为复杂的柱状图,并在图中进行适当标注

示例代码和结果如下:小疑问:如何自动选择合适的坐标位置对柱状图的数值大小进行标注?😂 clear; close all;x = 1:3;aa=[28.6321521955954 26.2453660695847 21.69102348512086.93747104431360 6.25442246899816 3.342835958564245.51365061796319 4.87

C# double[] 和Matlab数组MWArray[]转换

C# double[] 转换成MWArray[], 直接赋值就行             MWNumericArray[] ma = new MWNumericArray[4];             double[] dT = new double[] { 0 };             double[] dT1 = new double[] { 0,2 };

libsvm在matlab中的使用方法

原文地址:libsvm在matlab中的使用方法 作者: lwenqu_8lbsk 前段时间,gyp326曾在论坛里问libsvm如何在matlab中使用,我还奇怪,认为libsvm是C的程序,应该不能。没想到今天又有人问道,难道matlab真的能运行libsvm。我到官方网站看了下,原来,真的提供了matlab的使用接口。 接口下载在: http://www.csie.ntu.edu.

Matlab/Simulink中PMSM模型的反电动势系数和转矩系数

Matlab/Simulink中PMSM模型的反电动势系数和转矩系数_matlab pmsm-CSDN博客