快速学会有限元编程

《快速学会有限元编程》

一、前言

相信很多做过有限差分之后又想做做有限元的初学者会有和我一样的困惑，能看懂有限元算法的理论分析，但是真正应用到实际编程当中之前心里发怵，请教学过有限元程序的同学的时候，他们往往会，这个怎么怎么的简单，这个你怎么能不会？这个不就是什么什么吗bulabula...这时候你的心里一定和我一样有一万匹草泥马在心里奔腾，我特么的要是会还来问你，废话不多说，求人不如求己，这篇文章将会让你迅速掌握有限元最基础的编程思想。

二、以经典扩散方程为例

考虑如下扩散方程初边值问题

\begin{equation*}\label{eq1.1}
\left \{
\begin{array}{rl}
\frac{\partial u(x,t)}{\partial t}=\frac{\partial^2u(x,t)}{\partial t}+f(x,t), \quad & a< x< b,0< t\leq T, \\
u(a,t)=0,\quad u(b,t)=0,\quad & 0\leq t\leq T, \\
u(x,0)=\varphi (x),\quad & a\leq x \leq b, \\
\end{array}
\right.\tag{1}
\end{equation*}

首先我们需要对时间与空间区间进行剖分，其中$M,N$分别记为空间与时间方向的剖分，$h=\frac{b-a}{M}$, $ \tau = \frac{T}{N} $，$h,\tau$分别为空间与时间方向的步长，则$a=x_0<x_1<\cdots < x_M=b$, $0=t_0<t_1<\cdots< t_N=T$。

空间方向利用有限元逼近，可得（1）式的变分形式（弱形式）

$$\big( \frac{\partial u_h(t)}{\partial t},v_h \big)+a\big(u_h,v_h\big)=\big(f,v_h\big),\qquad v_n \in X_h \tag{2}$$.

其中时间方向利用向后差分离散,可得求解问题（1）的全离散格式

$$ \big( u^n_h,v_h \big)+\tau a\big(u^n_h,v_h\big)=\big( u^{n-1}_h,v_h \big)+\tau\big(f^n,v_h\big) \tag{3}$$

其中, 记双线形式$a(u^n_h,v_h)=\int_\Omega (u^n_h)‘(v_h)‘dx$为刚度矩阵（stiffness matrix），$(u^n_h,v_h)=\int_\Omega u^n_h v_h dx$为质量矩阵（mass matrix），$(f^n,v_h)=\int_\Omega f^n v_hdx$ 为载荷向量(load vector)。$\Omega=[a,b]=\underset{i=1,2,3,\cdots M-1}{\cup} \Omega_i$，$X_h=\{ \phi_i \}$为有限元空间，这里我们选用的是线性插值

$$\phi_i=\left \{
\begin{array}{ll}
\frac{x_{i+1}-x}{h},\qquad x\in \Omega_i,\\
\frac{x-x_{i-1}}{h},\qquad if\quad x\in \Omega_{i+1},\\
0, \qquad otherwise.\\
\end{array}
\right.\tag{4}$$

因此，我们需要获得的问题（1）的有限元解$u^n_h$可表示为

$$u^n_h=\sum\limits_{j=1}^{M-1} u_j^n \phi_j ,\tag{5}$$

其中$u^n_i$ 为问题（1）在点$(x_i,t_n)$处的数值解,将（5）式代入（3）式，我们可以获得求解问题（1）的有限元格式的代数方程组。.

$$\Big( \tau A + B \Big)u^n=Bu^{n-1}+\tau \widehat{f}.\tag{6}$$

其中$u^n=[u_1^n,u_2^n,\cdots,u_{M-1}^n]^T$, $A=(a_{ij})_{M-1\times M-1}$为刚度矩阵，由（4）式我们可以发现，对任意$i$, $\phi_i$只在$\Omega_i=[x_{i-1},x_i]$与$\Omega_{i+1}=[x_i,x_{i+1}]$两个区间不为零，其余区间即$|i-j|\geq 2$时$a_{ij}$均为零，则通过计算，我们发现

$$\begin{array}{rl}a_{i,i}&=\int_\Omega (\phi_i)‘(\phi_i)‘dx=\sum\limits_{i=1}^{M-1} \int_{\Omega_i} (\phi_i)‘(\phi_i)‘dx\\&=\int_{x_{i-1}}^{x_i} (\phi_i)‘(\phi_i)‘dx+\int_{x_i}^{x_{i+1}} (\phi_i)‘(\phi_i)‘dx\\&=\int_0^1 \frac{1}{h}d\xi+\int_0^1 \frac{1}{h} d\xi\\&=\frac{2}{h},\qquad i=1,2,\cdots, M-1.\tag{7}\end{array}$$

$$\begin{array}{rl}a_{i-1,i}=a_{i,i-1}&=\int_\Omega (\phi_i)‘(\phi_{i-1})‘dx=\sum\limits_{i=1}^{M-1} \int_{\Omega_i} (\phi_i)‘(\phi_{i-1})‘dx\\&=\int_{x_{i-1}}^{x_{i-1}} (\phi_i)‘(\phi_i)‘dx\\&=\int_0^1 \frac{-1}{h}d\xi\\&=\frac{-1}{h},\qquad i=2,3,\cdots,M-1.\tag{8}\end{array}$$

可知刚度矩阵$A$为一三对角矩阵，同理可知质量矩阵$B=(b_{ij})_{M-1\times M-1}$也为三对角矩阵，其矩阵元素由如下计算获得

$$\begin{array}{rl}b_{i,i}&=\int_\Omega (\phi_i)^2dx=\sum\limits_{i=1}^{M-1} \int_{\Omega_i} (\phi_i)^2dx\\&=\int_{x_{i-1}}^{x_i} (\phi_i)^2dx+\int_{x_i}^{x_{i+1}} (\phi_i)^2dx\\&=h \int_0^1 \xi^2d\xi+h\int_0^1 (1-\xi)^2d\xi\\&=\frac{2h}{3},\qquad i=1,2,\cdots, M-1.\tag{9}\end{array}$$

$$\begin{array}{rl}b_{i-1,i}=b_{i,i-1}&=\int_\Omega (\phi_i)(\phi_{i-1})dx=\sum\limits_{i=1}^{M-1} \int_{\Omega_i} (\phi_i)(\phi_{i-1})dx\\&=\int_{x_{i-1}}^{x_{i-1}} (\phi_i)(\phi_i)dx\\&=h\int_0^1 \xi(1-\xi)d\xi \\&=\frac{h}{6},\qquad i=2,3,\cdots,M-1.\tag{10}\end{array}$$

到了这里，我们终于可以松一口气了，主要的东西我们已经准备好了，，，，等等，，，还有载荷向量$\widehat{f}=[(f^n,\phi_1),(f^n,\phi_2),\cdots,(f^n,\phi_{M-1})]^T$还没计算，接下来给出载荷向量的计算

$$\begin{array}\Big( f^n,\phi_i \Big)&=\int_\Omega f^n,\phi_i dx=\sum\limits^{M-1}_{i=1} \int_{\Omega_{i}}f^n,\phi_i dx\\&=\int_{x_{i-1}}^{x_i} f(x)\phi_idx+\int_{x_i}^{x_{i+1}} f(x)\phi_idx\\&=h\int_0^1 f(x_{i-1}+h\xi)xdx+h\int_0^1 f(x_i+h\xi)(1-\xi)d\xi.\tag{11}\end{array}$$

为计算以上的一重积分，我在程序中采用的是17点$Gauss$求积公式，绝对符合我们对算法精度的要求！不会$Gauss$积分公式的同学可以查阅文献[2]pp116页.具体程序实现如下

function u = single_quad( a,b,f )
%   This rountine is written for Stiff/mass matrix terms, in which five-point
%   Gauss numerical integral is used.
%   f(x) is standar for integrand in interval [a,b]
%   g(x) is transformed from f(x), and [a,b] turn into [-1,1].
M = 17;
g = gauss_transform( f,a,b );
y(1)=0.1784841814958478558506775;y(2)=0.3512317634538763152971855;y(3)=0.5126905370864769678862466;
y(4)=0.6576711592166907658503022;y(5)=0.7815140038968014069252301;y(6)=0.8802391537269859021229557;
y(7)=0.9506755217687677612227170;y(8)=0.9905754753144173356754340;y(9)=0;y(10)=-y(8);y(11)=-y(7);
y(12)=-y(6);y(13)=-y(5);y(14)=-y(4);y(15)=-y(3);y(16)=-y(2);y(17)=-y(1);
A(1)=0.1765627053669926463252710;A(2)=0.1680041021564500445099707;A(3)= 0.1540457610768102880814316;
A(4)=0.1351363684685254732863200;A(5)=0.1118838471934039710947884;A(6)=0.0850361483171791808835354;
A(7)=0.0554595293739872011294402;A(8)=0.0241483028685479319601100;A(9)= 0.1794464703562065254582656;
A(10)=A(8);A(11)=A(7);A(12)=A(6);A(13)=A(5);A(14)=A(4);A(15)=A(3);A(16)=A(2);A(17)=A(1);
u = 0;

for i = 1:M
    u = u + A(i) * g( y(i) );
end

function g = gauss_transform( f,a,b )
        g = @(y) ( b - a ) / 2 * f( ( a+b ) / 2 + ( b-a ) / 2 * y );
end

end

到了这里，我们基本清楚了有限元格式基本算法的实现,以下我逐一介绍代码的实现过程。

在编程之前，我们自己需要构思一下，我们到底需要干些什么事情，第一步：输入算例的基本参数，第二步：区间剖分，第三步全离散格式（5）的封装，第四步：数据处理（收敛阶及误差），可视化。要想比较有条理的一步一步实现，我建议将2-4步单独封装成一个函数文件，第一步写成脚本文件，因为第一步需要我们手动输入模型参数，当你换另一个算例时只需要在第一步的文件中修改参数，接下来看看第一步:(取真解$u(x,t)=x^2(x-2)e^t$, $f(x,t)=(x^3-2x^2-6x+4)e^t$, $[a,b]=[0,2], T=1$, 初值$u(x,0)=x^2(x-2)$.)

脚本文件：run_main.m

clear;clc
% author: Dongdong Hu
% date: 16-Sep-2018
% version:8.3.0.532 (R2014a)
M = 100; N =100; Lx = 0; Rx = 2; Lt = 0; Rt = 1;
left_condation = @( t ) 0;
right_condation = @( t ) 0;
initial_condation = @( x ) x.^2 .* ( x-2 );
f = @( x,t ) exp( t ) .* ( x.^3 - 2 * x.^2 - 6 * x + 4 );
exact = @( x,t ) exp( t ) .* x.^2 .* ( x-2 );
show = show_solution(  );
LFEM = model_data( Lx,Rx,Lt,Rt,left_condation,right_condation,...
    initial_condation,exact);
[ x,t,u ] = Linear_FEM( M,N,LFEM,f );
ue = LFEM.exact( x,t );
show.image_3D( x,t,u );
show.image_2D( x,u,ue )
show.image_2D_error( x,abs( u-ue ) )
max_error = show.max_error( [ 10,20,40,80 ],[ 10,20,40,80 ].^2,LFEM,f );
max_rate = show.error_rate( max_error )
[ L2_error,space_step ] = show.L2_error( [ 10,20,40,80 ],[ 10,20,40,80 ].^2,LFEM,f );
L2_rate = show.error_rate( L2_error )
show.plot_rate( max_error,space_step );

　　第二步，我们需要进行网格剖分及初边值条件等已知参数的输入，这些都属于模型的固有属性

函数文件：model_data.m

function LFEM = model_data( Lx,Rx,Lt,Rt,left_condation,right_condation,...
    initial_condation,EXACT)
%MODEL_DATA 此处显示有关此函数的摘要
%   此处显示详细说明
LFEM = struct(‘space_grid‘,@space_grid,...
    ‘time_grid‘,@time_grid,‘left_boundary‘,@left_boundary,...
    ‘right_boundary‘,@right_boundary,‘initial_value‘,...
    @initial_value,‘exact‘,@exact);

    function [ x,h ] = space_grid( M )
            h = ( Rx - Lx ) / M;
            x = linspace( Lx,Rx,M+1 );
    end

    function [t,tau] = time_grid( N )
        tau = ( Rt - Lt ) / N;
        t = linspace( Lt,Rt,N+1 );
    end

    function u = left_boundary( t )
            u = left_condation( t );
    end

    function u = right_boundary( t )
            u = right_condation( t );
    end

    function u = initial_value( x )
            u = initial_condation( x );
    end

    function u = exact( x,t )
        u = zeros( length( x ),length( t ) );
        for k = 1:length( t )
                u( 1:end,k ) = EXACT( x,t( k ) );
        end
    end

end

　　第三步，我们需要将算法（5）封装

函数文件：Linear_FEM.m

function [ x,t,u ] = Linear_FEM( M,N,LFEM,f )
%L1_LINEAR_FEM Galerkin Finite elemeth method
%   this programming apply linear interpolation in finite elment method
%   with respect to x, the integer-order time partial derivative is
%   approximated by backward difference method.
u = zeros( M+1,N+1 );
[ x,h ] = LFEM.space_grid( M );
[t,tau] = LFEM.time_grid( N );

u( 1,1:N+1 ) = LFEM.left_boundary( t );
u( M+1,1:N+1 ) = LFEM.right_boundary( t );
u( 1:M+1,1 ) = LFEM.initial_value( x );

stiffness_matrix = zeros( M-1 );
mass_matrix = zeros( M-1 );

    for i = 1:M-1
        stiffness_matrix( i,i ) = 2 / h;
        mass_matrix( i,i ) = 2 / 3 * h;
        if i>1
            stiffness_matrix(i-1,i) = -1 / h;
            mass_matrix(i-1,i) = h / 6;
            stiffness_matrix(i,i-1) = -1 / h;
            mass_matrix(i,i-1) = h / 6;
        end
    end

    for n = 1:N
        nn = n + 1;

        load_vector = tau * h * single_quad( 0,1,@( xi ) f( x( 1:end-2 ) + h * xi,t( nn ) ) .* xi )‘...
            + tau * h * single_quad( 0,1,@( xi ) f( x( 2:end-1 ) + h * xi,t( nn ) ) .* ( 1-xi ) )‘...
            + mass_matrix * u( 2:end-1,nn-1 );
        u( 2:end-1,nn ) = ( tau * stiffness_matrix + mass_matrix ) \ ( load_vector );
    end

end

　　第四步，我们需要对我们所获得的数值解进行数据分析

函数文件：show_solution.m

function show = show_solution(  )
%SHOW_SOLUTION 此处显示有关此函数的摘要
%   此处显示详细说明
show = struct( ‘image_3D‘,@image_3D,‘image_2D‘,@image_2D,‘image_2D_error‘,@image_2D_error,...
    ‘max_error‘,@max_error,‘error_rate‘,@error_rate,‘L2_error‘,@L2_error,‘L2_rate‘,@L2_rate,...
    ‘plot_rate‘,@plot_rate );

    function image_3D( X,T,u )
       figure
       [ x,t ] = meshgrid( T,X );
       mesh( x,t,u );
       xlabel(‘x‘);ylabel(‘y‘);zlabel(‘u‘);
    end

    function image_2D( x,u,ue )
       figure
       plot( x,u(1:end,end),‘b*‘ );hold on;
       plot( x,ue( 1:end,end ),‘r-‘ );
       legend(‘numerical.sol‘,‘exact.sol‘);
       xlabel(‘x‘);ylabel(‘u( x,t=1 )‘);
    end

    function image_2D_error( x,u )
       figure
       plot( x,u(1:end,end) );
       xlabel(‘x‘);ylabel(‘|error( x,t=1 )|‘);
    end

    function Max_error = max_error( M,N,LFEM,f )
        Max_error = zeros( size( M ) );
        for i = 1:length( M )
           [ x,t,u ] = Linear_FEM( M(i),N(i),LFEM,f );
           ue = LFEM.exact( x,t );
           Max_error( i ) = max( max( abs( u-ue ) ) );
        end
    end

    function Error_rate = error_rate( max_error )
        Error_rate = log2( max_error(1:end-1) ./ max_error( 2:end ) );
    end

    function [ L2_Error,space_step ] = L2_error( M,N,LFEM,f )
        L2_Error = zeros( size( M ) );
        space_step = L2_Error;
        for i = 1:length( M )
          [ x,t,u ] = Linear_FEM( M(i),N(i),LFEM,f );
          ue = LFEM.exact( x,t );
          [ ~,h ] = LFEM.space_grid( M(i) );
          [~,tau] = LFEM.time_grid( N(i) );
          space_step( i ) = h;
          L2_Error(i) = sqrt( h * tau * sum( sum( ( u( 2:end-1, 2:end )-ue( 2:end-1, 2:end ) ).^2 ) ) );
        end
    end

    function plot_rate( error1,error2,space_step )
        figure
        loglog( space_step,error1,‘-kp‘ );hold on;
        loglog( space_step,error2,‘-r^‘ );hold on;
        loglog( space_step,space_step.^2,‘-g*‘ );
        legend(‘||error||_{\infty}‘,‘|| error ||_{L_2}‘,‘Ch^2‘);
        xlabel(‘log(1/h)‘); ylabel(‘error‘);
    end

end

　　到这里，整个编程到此结束，下面我来展示程序运行结果

原文地址：https://www.cnblogs.com/xtu-hudongdong/p/9657021.html

时间： 2025-01-30 01:51:42

快速学会有限元编程

快速学会有限元编程的相关文章

关键20小时，快速学会任何技能

怎样快速学会一门技术（转载）

为什么好多人想学Python 怎么快速学会高端技术

5个步骤快速学会自己建立个人网站

快速学会Python（实战）

Android 工程师如何快速学会web前段

快速学会Spring动态代理原理

[快速学会Swift第三方库] Eureka篇

5天学会jaxws-webservice编程第一天