萌新学习有限元方法解偏微分方程编程实例（三）

已于 2024-08-09 00:27:02 修改
阅读量980 收藏 12
点赞数 10
文章标签： 学习 matlab 数学建模 算法 线性代数 矩阵
于 2024-08-04 23:52:29 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_58189708/article/details/140913634
版权

一、题目

本人新手入门，从最简单最基础的例子学习，不足之处请大家多多批评指正。

以何晓明老师的课件例题为例：例题三和例题四

【本篇】博客只讨论一次线性的基函数下，三种边界条件和三种范数误差。

【下篇】博客将讨论二次多项式基函数下，三种边界条件和三种范数误差。

二、程序功能

1.主文件：run_1D_linear.m

%% 初始化
clc,clear
Bound = 1;  %选择边界条件，狄利克雷为1，纽曼为2，罗宾为3
N=32;       % N是有限元个数
xx=[0:1:N]; 
h=1.0/N;    % h是步长
xxx=xx.*h;  % xxx是刻度
nn=10;      % nn是积分精度
A=zeros(N+1,N+1);
B=zeros(N+1,1);

%% 古典解及其导数
uu  =  @(x) x.*cos(x);
duu =  @(x) cos(x) - x.*sin (x);

%% 逐个循环，生成单元刚度矩阵，并且组装总刚度矩阵
for j = 1:N
   x1=xxx(j);x2=xxx(j+1); 
   [AA,BB] = unit_1D_linear(N,h,j,x1,x2);
   A(j:j+1,j:j+1)=AA+A(j:j+1,j:j+1);
   B(j:j+1,  1  )=BB+B(j:j+1,  1  );
end

%% 边界条件 
switch fix(Bound)
    case 1
       [A,B] = Dirichlet_1D(A,B,uu);
    case 2
       [A,B] = Neu_1D(A,B,uu);
    case 3
       [A,B] = Robin_1D(A,B,uu);
end
  
%% 求解,作图
u = inv(A) * B;
uu(xxx);
plot(xxx,u,'b-*',xxx,uu(xxx),'r-')
legend('数值解(蓝色虚线)','真实解(红色实线)','Location', 'best');

%% 误差,L_2范数,L_wq,H1_semi范数,最大节点误差
[L_2,L_wq,H1_semi,max_err_node]= error_1D_linear(u,uu,duu,xxx,nn);

2.单元刚度/载荷矩阵：unit_1D_linear.m

function [AA,BB] = unit_1D_linear(N,h,j,x1,x2)

% unit_1D_linear.m生成单个单元刚度/载荷矩阵
% 请修改为自己的c(x) 
% AA是单元刚度矩阵
% BB是单元载荷矩阵

c =  @(x)   exp(x)*(1/h)*(1/h);
FB1 = @(x) -exp(x).*(cos(x)-2.*sin(x)-x.*cos(x)-x.*sin(x)).*(x2-x)*(1/h);
FB2 = @(x) -exp(x).*(cos(x)-2.*sin(x)-x.*cos(x)-x.*sin(x)).*(x-x1)*(1/h); 
 
AA=zeros(2,2);
BB=zeros(2,1);
AA(1,1)=  integral(c,x1,x2);
AA(1,2)= -integral(c,x1,x2);
AA(2,1)= -integral(c,x1,x2);
AA(2,2)=  integral(c,x1,x2);

BB(1,1)=  integral(FB1,x1,x2);
BB(2,1)=  integral(FB2,x1,x2);

end

3.边界条件的处理

（1）Dirichlet_1D.m  处理狄利克雷边界条件

function [A,B] = Dirichlet_1D(A,B,uu) 
   % 一维问题的狄利克雷边界条件
   A(1,:)=0;
   A(1,1)=1;
   A(end,:)=0;
   A(end,end)=1; 
   B(1)=uu(0);
   B(end)=uu(1);
end

（2）Neu_1D.m  处理纽曼边界条件

function [A,B] = Neu_1D(A,B,uu)
% 一维问题的纽曼边界条件
   A(1,:)=0;
   A(1,1)=1;
   B(1)=uu(0);
   %注意，因为我们选择一次线性“帽子”函数
   %只有z最后一个基函数的内积不为零，在b点值为1，而不是为0
   B(end)=B(end)+(cos(1)-sin(1)).*exp(1).*1  ;
end

（3）Robin_1D.m 处理罗宾边界条件

function [A,B] = Robin_1D(A,B,uu)
  %一维问题的罗宾边界条件
   A(end,:)=0; 
   A(end,end)=1; 
   B(end)=uu(1);  
   %注意，因为我们选择一次线性“帽子”函数
   %只有第一个基函数及其内积，在a点值为1，而不是为0
   A(1,1) = A(1,1)+1.*exp(0).*1.*1;
   B(1)  =  B(1)-1.*exp(0).*1; 
end

4.误差分析：error_1D_linear.m

function [L_2,L_wq,H1_semi,max_err_node]= error_1D_linear(u,uu,duu,Mesh,nn)
%% 节点误差
max_err_node = max(abs(u-uu(Mesh)'));%无穷范数为例

%% L_2范数,L_wq,H1_semi范数
L_2  = 0;
H1_semi =0;
L_wq = 0;
L_wq_Max = [2,1:size(Mesh,2)-1];

for k = 1 : size(Mesh,2)-1
    % 每个有限元上，只有两个“帽子”函数(基函数)不为零
    % L_2范数,用Uh乘以基函数，使数值解连续
    % L_wq范数找每个区间最大差值，再找整个定义域上的最大值
    u_L2       =  @(x) -u( k ).*( x - Mesh(k+1) )./(Mesh(k+1) - Mesh(k))  ...
                       +u(k+1).*( x - Mesh( k ) )./(Mesh(k+1) - Mesh(k));  
    %注意，fminbn取的反而是极小值，绝对值加了负号
    uu_L       =  @(x) -abs ( feval(u_L2,x) - feval(uu,x) ) ;
    uuu_L2     =  @(x)      ( feval(u_L2,x) - feval(uu,x) ).^2;
    u_Hhalf    =  @(x)  (u(k+1) - u( k ))./(Mesh(k+1) - Mesh(k));
    % H1_semi范数同理
    uu_Hhalf   =  @(x)  ( feval(u_Hhalf,x) - feval(duu,x) ).^2;
    
    %存储每个区间的极大值
    %注意，fminbn取的反而是极小值，绝对值加了负号
    [L_wq_Max(1,k) , L_wq_Max(2,k)] = fminbnd(uu_L,Mesh(k),Mesh(k+1));
    %分段数值积分
    L_2     = L_2    + mysimp( uuu_L2   , Mesh(k+1) , Mesh(k) , nn);
    H1_semi = H1_semi + mysimp( uu_Hhalf , Mesh(k+1) , Mesh(k) , nn);
end 

  %L_2和H1_semi开根号
  L_2     = sqrt( L_2 );
  H1_semi = sqrt(H1_semi);
  %L_wq取最大值，注意，fminbn取的反而是极小值，这里去掉负号
  L_wq = max  ( -L_wq_Max(2,:));
  
  %% 展示输出
  format long;
  fprintf('L_2范数：');
  L_2
  fprintf('L_wq范数：');
  L_wq
  fprintf('H1_semi范数：');
  H1_semi
  fprintf('节点最大误差：')
  max_err_node
  
end

5.辛普森积分：mysimp.m

function I = mysimp(f,b,a,nn)
% 注释
% nn是子区间的个数,可以调节的精度；
%f为被积函数；
%a,b分别为积分的上下限；
%I是梯形总面积，即所求积分数值；
%feval('f',x)是泛函,例如sin(x),
 
%% 辛普森公式求积分
h=(b-a)/(2*nn);
s1=0;
s2=0;
for k=1:nn
    x=a+h*(2*k-1);
    s1=s1+feval(f,x);
end
for k=1:(nn-1)
    x=a+h*2*k;
    s2=s2+feval(f,x);
end
I = h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3;
 
end

三、结果与分析

1.对比三种边界条件的误差(h=4)

（1）狄利克雷边界条件

（2）纽曼边界条件

（3）罗宾边界条件

2.范数误差

（1）双对数图

（2）误差数值

3.小结

（1）我的L2和H1半范的结果都与标准答案精确符合，我的L∞范数与标准答案保持在同一个数量级，但是略大10%左右，通过图像分析，我猜想可能是我使用MATLAB内置函数fminbnd寻找极值，找到的最大误差可能更大，导致出现了一些偏差，但不清楚具体原因。

（2）为什么改变边界条件后，L2和L∞范数具有较大改变，H1半范几乎不变？请高手指点。

（3）有没有更简易的程序格式，更方便地调用基函数或导数？请高手指点？

shx48716
关注
偏微分方程的有限元解法及可视化（使用VTK）
FdmPatch的博客
09-22 373
有限元法的基本思想是将求解域划分为离散的小单元，然后在每个单元内构建近似解。对于每个单元，我们引入一个试探函数（Test Function）和一个权重函数（Weight Function），并将泊松方程转化为弱形式。通过对所有单元的求和，我们可以得到离散的系统方程。在本文中，我们选择经典的泊松方程（Poisson Equation）作为示例。其中，u是未知函数，表示待求解的物理量，f是已知函数，表示外部源项。接下来，我们将使用Python代码实现有限元法求解泊松方程，并使用VTK进行可视化。
二阶抛物型偏微分方程有限元法matlab,matlab使用有限元方法求解偏微分方程
weixin_34929072的博客
03-17 1465
【实例简介】大量的有限元法求解偏微分方程的实例程序，程序里有详细的解释语句【实例截图】【核心代码】fem_matlab大量的有限元法求解偏微分方程的实例程序(程序有解释)└── fem_matlab大量的有限元法求解偏微分方程的实例程序(程序有解释)├── fem_matlab│ └── fem_matlab│ ├── Ex10101.m│ ├── Ex10102.m│...
偏微分方程的有限元解法
feyily的博客
04-17 4551
本文我们以求解泊松方程为例，讲述微分方程的有限元解法，附FEniCS代码 方程定义 首先是泊松方程的表达式 −Δu=f in Ω(1) -\Delta u=f \qquad \text{ in } \Omega \tag{1} −Δu=f in Ω(1) 其边界条件 u∣ΓD=u0 on ΓD∂u∂n∣ΓN=g on ΓN=∂Ω−ΓD u|_{\Gamma_D}=u_0 \qquad \text{ on } \Gamma_D.
萌新学习有限元方法解偏微分方程编程实例（一）
weixin_58189708的博客
07-29 291
萌新学习有限元方法解偏微分方程编程实例（一）
偏微分方程定解问题的有限元方法（Matlab）
08-30
有限元方法解PDE，二维边值问题，矩形剖分。从区域剖分、计算刚度矩阵，到解转化的代数方程，比较详细！
matlab使用有限元方法求解偏微分方程
04-08
matlab使用有限元方法求解偏微分方程
萌新学习有限元方法解偏微分方程编程实例（二）
最新发布
weixin_58189708的博客
08-01 468
本人新手入门，尝试从最简单最基础的例子编写程序，利用二次多项式函数作为基底，求解带有狄利克雷边界问题的一维偏微分方程。
matlab使用有限元方法求解偏微分方程.zip
07-29
matlab使用有限元方法求解偏微分方程
Matlab-PDE工具箱有限元法求解偏微分方程
05-07
此文档用于了解Matlab-PDE工具箱有限元法求解偏微分方程的具体步骤
有限元法求解常微分方程--两点边值问题
08-22
有限元法求解常微分方程 方程类型为 -u''(x)+q*u=f(x),u(a)=0,u'(b)=0,x \in (a,b) 其中q为常数 数值分析程序
有限元法基础
04-24
对于学习有限元原理有很大帮助，能够刚初学者更好的理解
求解偏微分方程开源有限元软件deal.II学习--Step 48
摄心神，致虚极，守静笃
01-22 2007
求解偏微分方程开源有限元软件deal.II学习--Step 48 Posted on 2016-11-29   |   In computational material science   |   暂无评论 引子 本例提供了一个框架来应用MatrixFree类，既包括求解非线性偏微分方程过程，同时演示MatrixFree类怎样处理“constraints”以及如何在分布式
有限元入门： 用变分法求解一维常微分方程(二)
waitingwinter的博客
05-20 2816
承接上一篇博客 https://blog.csdn.net/waitingwinter/article/details/106164350 为了使结构完整，我们再次给出所求解的问题: {−U′′(t)+U(t)=(π2+1)sin⁡πt,U(0)=U(1)=0,\left\{ \begin{aligned} & -U''(t) +U(t) = (\pi^2 + 1)\sin \pi t,\\ & U(0)=U(1) =0, \end{aligned} \right.{​−U′′(t)+U(t
有限元微分方程求解方法，能量原理，瑞利里兹法，伽辽金法（曾攀有限元分析）
热门推荐
瓜地考拉
04-16 1万+
1微分方程求解方法 以一个承受水平均布载荷的拉杆，计算其位移场u(x)为例，介绍有限元微分方程常用的求解方法。\qquad以一个承受水平均布载荷的拉杆，计算其位移场u_(x)为例，\\ 介绍有限元微分方程常用的求解方法。以一个承受水平均布载荷的拉杆，计算其位移场u(​x)为例，介绍有限元微分方程常用的求解方法。 1.1 有限差分法 微分方程中的高阶项用高阶差商代替，高阶差商使用低阶差商层层嵌套表...
有限元法求解微分方程
11-21
有限元法是一种数值解微分方程的方法，它将求解区域分割成许多小的单元，每个单元内部的微分方程可以用简单的代数方程来近似表示。通过将所有单元的代数方程组合起来，可以得到整个求解区域的微分方程的数值解。以下是有限元法求解微分方程的一般步骤： 1. 将求解区域分割成许多小的单元，每个单元内部的微分方程可以用简单的代数方程来近似表示。 2. 对于每个单元，选择一个适当的插值函数来近似表示未知函数。 3. 将每个单元的代数方程组合起来，得到整个求解区域的微分方程的数值解。 4. 对于边界条件，可以采用类似插值函数的方法来近似表示。 5. 通过求解代数方程组得到数值解。 以下是一个使用有限元法求解二阶常微分方程的Python代码示例： ```python from scipy.sparse import lil_matrix from scipy.sparse.linalg import spsolve import numpy as np def fem_solve(a, b, c, f, xa, xb, n): h = (xb - xa) / n x = np.linspace(xa, xb, n+1) A = lil_matrix((n+1, n+1)) b = np.zeros(n+1) for i in range(1, n): A[i, i-1] = a(x[i]) / h**2 - b(x[i]) / (2*h) A[i, i] = -2 * a(x[i]) / h**2 + c(x[i]) A[i, i+1] = a(x[i]) / h**2 + b(x[i]) / (2*h) b[i] = f(x[i]) A[0, 0] = 1 A[n, n] = 1 b[0] = xa b[n] = xb u = spsolve(A.tocsr(), b) return x, u # 示例：求解 y'' + y = x, y(0) = 0, y(1) = 1 a = lambda x: 1 b = lambda x: 0 c = lambda x: 1 f = lambda x: x xa = 0 xb = 1 n = 10 x, u = fem_solve(a, b, c, f, xa, xb, n) print(u) ```
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值