两点边值问题的有限元方法以及边值条件的处理示例

news2024/11/23 15:11:03

文章目录

  • 引言
    • 题目
    • 补全方程
    • 刚度矩阵
    • 构造基底
    • 边值条件非齐次
      • 左边值条件非齐次
      • 右边值条件非齐次
      • 非齐次边值条件有限元方程
    • 求数值解
      • 直接求总刚度矩阵
      • 先求单元刚度矩阵

引言

本文参考李荣华教授的《偏微分方程数值解法》一书

题目

对于非齐次第二边值问题
{ − d d x ( p d u d x ) + q u = f , a < x < b u ( a ) = α , u ′ ( b ) = β \begin{cases} -\dfrac{d}{dx}(p\dfrac{du}{dx})+qu=f,a<x<b\\ u(a)=\alpha,u'(b)=\beta \end{cases} dxd(pdxdu)+qu=f,a<x<bu(a)=α,u(b)=β
已知精确解 u ( x ) = sin ⁡ 5 x + x 3 ( 2 − x ) + 2 u(x)=\sin5x+x^3(2-x)+2 u(x)=sin5x+x3(2x)+2,并且
{ p ( x ) = cos ⁡ x q ( x ) = x \begin{cases} p(x)=\cos x\\ q(x)=x \end{cases} {p(x)=cosxq(x)=x
求剖分区间数分别为16和32时的精确解并画图

补全方程

使用 S a g e M a t h SageMath SageMath进行符号运算,求出 u ′ ( x ) u'(x) u(x) f ( x ) f(x) f(x),代码如下

x = var('x')
u = sin(5*x)+x^3*(2-x)+2
p = cos(x)
q = x
f = -derivative(p * derivative(u, x), x) + q * u
print(f)
print(derivative(u, x))

结果为:

-((x - 2)*x^3 - sin(5*x) - 2)*x + (6*(x - 2)*x + 6*x^2 + 25*sin(5*x))*cos(x) - (3*(x - 2)*x^2 + x^3 - 5*cos(5*x))*sin(x)
-3*(x - 2)*x^2 - x^3 + 5*cos(5*x)

于是原方程为
{ − d d x ( p d u d x ) + q u = − ( ( x − 2 ) x 3 − sin ⁡ 5 x − 2 ) x + ( 6 ( x − 2 ) x + 6 x 2 + 25 sin ⁡ 5 x ) cos ⁡ ( x ) − ( 3 ( x − 2 ) x 2 + x 3 − 5 cos ⁡ 5 x ) sin ⁡ x , a < x < b u ( a ) = sin ⁡ 5 a + a 3 ( 2 − a ) + 2 u ′ ( b ) = − 3 ( x − 2 ) x 2 − x 3 + 5 cos ⁡ 5 x \begin{cases} -\dfrac{d}{dx}(p\dfrac{du}{dx})+qu=-((x - 2)x^3 - \sin5x - 2)x +\\ (6(x - 2)x + 6x^2 + 25\sin 5x)\cos(x) - (3(x - 2)x^2 + x^3 - 5\cos 5x)\sin x,\\&a<x<b\\ u(a)=\sin5a+a^3(2-a)+2\\ u'(b)=-3(x - 2)x^2 - x^3 + 5\cos5x \end{cases} dxd(pdxdu)+qu=((x2)x3sin5x2)x+(6(x2)x+6x2+25sin5x)cos(x)(3(x2)x2+x35cos5x)sinx,u(a)=sin5a+a3(2a)+2u(b)=3(x2)x2x3+5cos5xa<x<b

刚度矩阵

对于网格剖分
a = x 0 < x 1 < ⋯ < x n = b a=x_0<x_1<\cdots<x_n=b a=x0<x1<<xn=b
对于每个单元 I i = [ x i − 1 , x i ] I_i = [x_{i-1}, x_i] Ii=[xi1,xi],按照线性插值公式,有数值解
u h = x i − x h i u i − 1 + x − x i − 1 h i u i , x ∈ I i , i = 1 , 2 , ⋯   , n u_h = \frac{x_i-x}{h_i}u_{i-1} + \frac{x-x_{i-1}}{h_i}u_i,x\in I_i,i=1,2,\cdots,n uh=hixixui1+hixxi1ui,xIi,i=1,2,,n
其中 h i = x i − x i − 1 h_i=x_i-x_{i-1} hi=xixi1

为使按段插值标准化,通常用仿射变换
ξ = x − x i − 1 h i \xi =\frac{x-x_{i-1}}{h_i} ξ=hixxi1
I i I_i Ii变到 ξ \xi ξ轴上的参考单元 [ 0 , 1 ] [0,1] [0,1],令
N 0 ( ξ ) = 1 − ξ , N i ( ξ ) = ξ N_0(\xi)=1-\xi,N_i(\xi)=\xi N0(ξ)=1ξ,Ni(ξ)=ξ

u h = N 0 ( ξ ) u i − 1 + N 1 ( ξ ) u i u_h=N_0(\xi)u_{i-1}+N_1(\xi)u_i uh=N0(ξ)ui1+N1(ξ)ui
代入泛函
J ( u ) = 1 2 ∫ a b ( p u ′ 2 + q u 2 − 2 f u ) d x J(u)=\frac{1}{2}\int_a^b (pu'^2+qu^2-2fu)dx J(u)=21ab(pu′2+qu22fu)dx

KaTeX parse error: Undefined control sequence: \part at position 8: \frac{\̲p̲a̲r̲t̲ ̲J(u_h)}{\part u…
其中, u j = u h ( x j ) , j = 1 , 2 , ⋯   , n , u 0 = 0 u_j=u_h(x_j),j=1,2,\cdots,n,u_0=0 uj=uh(xj),j=1,2,,n,u0=0
{ a j − 1 , j = ∫ 0 1 [ − p ( x j − 1 + h j ξ ) / h j + h j q ( x j − 1 + h j ξ ) ξ ( 1 − ξ ) ] d ξ a j + 1 , j = ∫ 0 1 [ − p ( x j + h j + 1 ξ ) / h j + 1 + h j + 1 q ( x j + h j + 1 ξ ) ξ ( 1 − ξ ) ] d ξ a j j = ∫ 0 1 [ p ( x j − 1 + h j ξ ) / h j + h j q ( x j − 1 + h j ξ ) ξ 2 ] d ξ + ∫ 0 1 [ p ( x j + h j + 1 ξ ) / h j + 1 + h j + 1 q ( x j + h j + 1 ξ ) ( 1 − ξ ) 2 ] d ξ b j = ∫ 0 1 [ h j f ( x j − 1 + h j ξ ) ξ ] d ξ + ∫ 0 1 [ h j + 1 f ( x j + h j + 1 ξ ) ( 1 − ξ ) ] d ξ ⋯ ( ∗ ) \begin{cases} a_{j-1,j}&=\int_0^1[-p(x_{j-1}+h_j\xi)/h_j+h_jq(x_{j-1}+h_j\xi)\xi(1-\xi)]d\xi\\ a_{j+1,j}&=\int_0^1[-p(x_{j}+h_{j+1}\xi)/h_{j+1}+h_{j+1}q(x_{j}+h_{j+1}\xi)\xi(1-\xi)]d\xi\\ a_{jj}&=\int_0^1[p(x_{j-1}+h_j\xi)/h_j+h_jq(x_{j-1}+h_j\xi)\xi^2]d\xi+\\ &\int_0^1[p(x_{j}+h_{j+1}\xi)/h_{j+1}+h_{j+1}q(x_{j}+h_{j+1}\xi)(1-\xi)^2]d\xi\\ b_{j}&=\int_0^1[h_jf(x_{j-1}+h_j\xi)\xi]d\xi+\int_0^1[h_{j+1}f(x_{j}+h_{j+1}\xi)(1-\xi)]d\xi \end{cases}\cdots(*) aj1,jaj+1,jajjbj=01[p(xj1+hjξ)/hj+hjq(xj1+hjξ)ξ(1ξ)]dξ=01[p(xj+hj+1ξ)/hj+1+hj+1q(xj+hj+1ξ)ξ(1ξ)]dξ=01[p(xj1+hjξ)/hj+hjq(xj1+hjξ)ξ2]dξ+01[p(xj+hj+1ξ)/hj+1+hj+1q(xj+hj+1ξ)(1ξ)2]dξ=01[hjf(xj1+hjξ)ξ]dξ+01[hj+1f(xj+hj+1ξ)(1ξ)]dξ()
于是可以得到总刚度矩阵 K K K,满足
K u = b Ku=b Ku=b
发现 a j j a_{jj} ajj b j b_j bj由不同的两项相加,于是可以得到单元刚度矩阵 K ( i ) K^{(i)} K(i)
K ( i ) = ( a i − 1 , i − 1 ( i ) a i − 1 , i ( i ) a i , i − 1 ( i ) a i i ( i ) ) K^{(i)}= \begin{pmatrix} a^{(i)}_{i-1,i-1}&a^{(i)}_{i-1,i}\\ a^{(i)}_{i,i-1}&a^{(i)}_{ii} \end{pmatrix} K(i)=(ai1,i1(i)ai,i1(i)ai1,i(i)aii(i))
其中
{ a i − 1 , i − 1 ( i ) = a i − 1 , i − 1 ( i ) = ∫ 0 1 [ − p ( x i − 1 + h i ξ ) / h i + h i q ( x i − 1 + h i ξ ) ξ ( 1 − ξ ) ] d ξ a i − 1 , i − 1 ( i ) = ∫ 0 1 [ p ( x i − 1 + h i ξ ) / h i + h i q ( x i − 1 + h i ξ ) ( 1 − ξ ) 2 ] d ξ a i i ( i ) = ∫ 0 1 [ p ( x i + h i + 1 ξ ) / h i + 1 + h i + 1 q ( x i + h i + 1 ξ ) ξ 2 ] d ξ ⋯ ( ∗ ∗ ) \begin{cases} a^{(i)}_{i-1,i-1}=a^{(i)}_{i-1,i-1}=\int_0^1[-p(x_{i-1}+h_i\xi)/h_i+h_iq(x_{i-1}+h_i\xi)\xi(1-\xi)]d\xi\\ a^{(i)}_{i-1,i-1}=\int_0^1[p(x_{i-1}+h_i\xi)/h_i+h_iq(x_{i-1}+h_i\xi)(1-\xi)^2]d\xi\\ a^{(i)}_{ii}=\int_0^1[p(x_{i}+h_{i+1}\xi)/h_{i+1}+h_{i+1}q(x_{i}+h_{i+1}\xi)\xi^2]d\xi\\ \end{cases}\cdots (**) ai1,i1(i)=ai1,i1(i)=01[p(xi1+hiξ)/hi+hiq(xi1+hiξ)ξ(1ξ)]dξai1,i1(i)=01[p(xi1+hiξ)/hi+hiq(xi1+hiξ)(1ξ)2]dξaii(i)=01[p(xi+hi+1ξ)/hi+1+hi+1q(xi+hi+1ξ)ξ2]dξ()
不难发现,总刚度矩阵可以由单元刚度矩阵 K ( i ) K^{(i)} K(i)“组装”而成,即对于刚度矩阵 K K K中的元素 a i j a_{ij} aij和单元刚度矩阵 K ( i ) K^{(i)} K(i)中的元素 a i j ( i ) a_{ij}^{(i)} aij(i),有如下关系
a i j = { a i , i − 1 ( i ) , j = i − 1 a i i ( i ) + a i i ( i + 1 ) , j = i a i , i + 1 ( i + 1 ) , j = i + 1 0 , ∣ j − i ∣ ≥ 2 a_{ij}= \begin{cases} a_{i,i-1}^{(i)}, &j=i-1\\ a_{ii}^{(i)}+a_{ii}^{(i+1)}, &j=i\\ a_{i,i+1}^{(i+1)}, &j=i+1\\ 0, &|j-i|\ge 2 \end{cases} aij= ai,i1(i),aii(i)+aii(i+1),ai,i+1(i+1),0,j=i1j=ij=i+1ji2
不难发现,总刚度矩阵 K K K为三对角矩阵

对应的,向量 b b b也可以进行“拆分”,令
{ f i − 1 ( i ) = ∫ 0 1 [ h j f ( x j − 1 + h j ξ ) ( 1 − ξ ) ] d ξ f i ( i ) = ∫ 0 1 [ h j + 1 f ( x j + h j + 1 ξ ) ξ ] d ξ \begin{cases} f^{(i)}_{i-1}=\int_0^1[h_jf(x_{j-1}+h_j\xi)(1-\xi)]d\xi\\ f^{(i)}_{i}=\int_0^1[h_{j+1}f(x_{j}+h_{j+1}\xi)\xi]d\xi \end{cases} {fi1(i)=01[hjf(xj1+hjξ)(1ξ)]dξfi(i)=01[hj+1f(xj+hj+1ξ)ξ]dξ
则有
b 1 = f 1 ( 1 ) + f 1 ( 2 ) , b 1 = f 2 ( 2 ) + f 2 ( 3 ) , ⋯   , b 1 = f n ( n ) ⋯ ( ∗ ∗ ∗ ) b_1=f_1^{(1)}+f_1^{(2)},b_1=f_2^{(2)}+f_2^{(3)},\cdots,b_1=f_n^{(n)}\cdots (***) b1=f1(1)+f1(2),b1=f2(2)+f2(3),,b1=fn(n)()

构造基底

再单元 I i , I i + 1 I_i,I_{i+1} Ii,Ii+1考察线性插值公式及 u i u_i ui的系数,我们自然对每一节点 x i x_i xi构造山形函数:
{ φ i ( x ) = { 1 + x − x i h i , x i − 1 ≤ x ≤ x i 1 − x − x i h i + 1 , x i ≤ x ≤ x i + 1 0 , o t h e r w i s e i = 1 , 2 , ⋯   , n − 1 φ n ( x ) = { 1 + x − x n h n , x n − 1 ≤ x ≤ x n 0 , o t h e r w i s e \begin{cases} \varphi_i(x)&=\begin{cases} 1+\dfrac{x-x_i}{h_i},&x_{i-1}\le x\le x_i\\ 1-\dfrac{x-x_i}{h_{i+1}},&x_{i}\le x\le x_{i+1}\\ 0,&otherwise \end{cases}\\ &i=1,2,\cdots,n-1\\ \varphi_n(x)&=\begin{cases} 1+\dfrac{x-x_n}{h_n},&x_{n-1}\le x\le x_n\\ 0,&otherwise \end{cases}\\ \end{cases} φi(x)φn(x)= 1+hixxi,1hi+1xxi,0,xi1xxixixxi+1otherwisei=1,2,,n1= 1+hnxxn,0,xn1xxnotherwise
于是 u h ( x ) = ∑ i = 1 n u i φ i ( x ) , u i = u h ( x i ) u_h(x)=\sum\limits_{i=1}^{n}u_i\varphi_i(x),u_i=u_h(x_i) uh(x)=i=1nuiφi(x),ui=uh(xi)

边值条件非齐次

左边值条件非齐次

u ( α ) = a u(\alpha)=a u(α)=a

则增加一基函数
φ 0 ( x ) = { a − x − x 0 h 1 , x 0 ≤ x ≤ x 1 0 , o t h e r w i s e \varphi_0(x)=\begin{cases} a-\dfrac{x-x_0}{h_1},&x_0\le x\le x_1\\ 0,&otherwise \end{cases} φ0(x)= ah1xx0,0,x0xx1otherwise
u h u_h uh表为
u h = α φ 0 ( x ) + ∑ i = 1 n u i φ i ( x ) u_h=\alpha \varphi_0(x)+\sum_{i=1}^n u_i\varphi_i(x) uh=αφ0(x)+i=1nuiφi(x)
有限元方程为
∑ i = 1 n a ( φ i , φ j ) u i = ( f , φ j ) − α a ( φ 0 , φ j ) , j = 1 , 2 , ⋯   , n \sum_{i=1}^n a(\varphi_i,\varphi_j)u_i=(f,\varphi_j)-\alpha a(\varphi_0,\varphi_j),j=1,2,\cdots,n i=1na(φi,φj)ui=(f,φj)αa(φ0,φj),j=1,2,,n
实际上只修改了第一个方程的右端,因为当 j = 2 , 3 , ⋯   , n j=2,3,\cdots,n j=2,3,,n时, a ( φ 0 , φ j ) = 0 a(\varphi_0,\varphi_j)=0 a(φ0,φj)=0

右边值条件非齐次

则右端修改为 ( f , φ j ) + p ( b ) β φ j ( b ) (f,\varphi_j)+p(b)\beta\varphi_j(b) (f,φj)+p(b)βφj(b),有限元方程变为
∑ i = 1 n u i ∫ a b [ p φ i ′ φ j ′ + q φ i φ j ] d x = ( f , φ j ) + p ( b ) β φ j ( b ) , j = 1 , 2 , ⋯   , n \sum_{i=1}^n u_i\int_a^b[p\varphi_i'\varphi_j'+q\varphi_i\varphi_j]dx=(f,\varphi_j)+p(b)\beta\varphi_j(b),j=1,2,\cdots,n i=1nuiab[pφiφj+qφiφj]dx=(f,φj)+p(b)βφj(b),j=1,2,,n
实际上只修改了第 n n n个方程的右端项,因为当 j = 1 , 2 , ⋯   , n − 1 j=1,2,\cdots,n-1 j=1,2,,n1时, φ j ( b ) = 0 \varphi_j(b)=0 φj(b)=0

非齐次边值条件有限元方程

综上,我们得到非齐次边值条件有限元方程
( a ( φ 1 , φ 1 ) ⋯ a ( φ n , φ 1 ) a ( φ 1 , φ 2 ) ⋯ a ( φ n , φ 2 ) ⋮ ⋱ ⋮ a ( φ 1 , φ n − 1 ) ⋯ a ( φ n , φ n − 1 ) a ( φ 1 , φ n ) ⋯ a ( φ n , φ n ) ) ( u 1 u 2 ⋮ u n − 1 u n ) = ( ( f , φ 1 ) − α a ( φ 0 , φ 1 ) ( f , φ 2 ) ⋮ ( f , φ n − 1 ) ( f , φ n ) + p ( b ) β φ n ( b ) ) \begin{pmatrix} a(\varphi_1,\varphi_1)&\cdots&a(\varphi_n,\varphi_1)\\ a(\varphi_1,\varphi_2)&\cdots&a(\varphi_n,\varphi_2)\\ \vdots&\ddots&\vdots\\ a(\varphi_1,\varphi_{n-1})&\cdots&a(\varphi_n,\varphi_{n-1})\\ a(\varphi_1,\varphi_n)&\cdots&a(\varphi_n,\varphi_n)\\ \end{pmatrix} \begin{pmatrix} u_1\\ u_2\\ \vdots\\ u_{n-1}\\ u_n \end{pmatrix} =\begin{pmatrix} (f,\varphi_1)-\alpha a(\varphi_0,\varphi_1)\\ (f,\varphi_2)\\ \vdots\\ (f,\varphi_{n-1})\\ (f,\varphi_n)+p(b)\beta\varphi_n(b) \end{pmatrix} a(φ1,φ1)a(φ1,φ2)a(φ1,φn1)a(φ1,φn)a(φn,φ1)a(φn,φ2)a(φn,φn1)a(φn,φn) u1u2un1un = (f,φ1)αa(φ0,φ1)(f,φ2)(f,φn1)(f,φn)+p(b)βφn(b)

求数值解

我们以区间 [ − 1 , 2 ] [-1,2] [1,2]为例做均匀剖分

在求解数值解 u = ( u 1 , u 2 , ⋯   , u n ) T u=(u_1,u_2,\cdots,u_n)^T u=(u1,u2,,un)T时,我们可以直接求总刚度矩阵,也可以先求单元刚度矩阵再进行“组装”

直接求总刚度矩阵

M a t l a b Matlab Matlab代码如下

% 直接求刚度矩阵
n = ;
a = -1;
b_ = 2;

A = zeros(n, n);
b = zeros(n, 1);

p = zeros(1, n+1);
I = zeros(2, n);
for i=1:1:n+1
    p(1,i) = a + (b_-a)/n.*(i-1);
    if i<n+1
        I(1,i) = i;
        I(2,i) = i+1;
    end
end

u0 = @(x) sin(5.*x)+x.^3.*(2-x)+2;
du = @(x) -3.*(x - 2).*x.^2 - x.^3 + 5.*cos(5.*x);
alpha = u0(a);
beta = du(b_);
px = @(x) cos(x);
qx = @(x) x;
f = @(x) -((x - 2).*x.^3 - sin(5*x) - 2).*x + (6*(x - 2).*x + 6*x.^2 + 25*sin(5*x)).*cos(x) - (3*(x - 2).*x.^2 + x.^3 - 5*cos(5*x)).*sin(x);
aii = @(x1, x2, h) integral(@(ksi) (px(x1+h.*ksi)/h + h.*qx(x1+h.*ksi).*ksi.^2 + px(x2+h.*ksi)/h + h.*qx(x2+h.*ksi).*(1-ksi).^2), 0, 1);
a23 = @(x, h) integral(@(ksi) (-px(x+h.*ksi)/h + h.*qx(x+h.*ksi).*ksi.*(1-ksi)),0, 1);
ann = @(x, h) integral(@(ksi) (px(x+h.*ksi)/h + h.*qx(x+h.*ksi).*ksi.^2), 0, 1);
bi = @(x1, x2, h) integral(@(ksi) (h.*f(x1+h.*ksi).*ksi + h.*f(x2+h.*ksi).*(1-ksi)), 0, 1);
bn = @(x, h) integral(@(ksi) (h.*f(x+h.*ksi).*ksi), 0, 1);

for i=1:1:n-1
    x1 = p(1, I(1,i));
    x2 = p(1, I(2,i));
    h = x2 - x1;
    A(I(1,i), I(1,i)) = aii(x1, x2, h);
    b(I(1, i), 1) = bi(x1, x2, h);
end
for i=1:1:n-1
    x1 = p(1, I(1,i+1));
    x2 = p(1, I(2,i+1));
    h = x2 - x1;
    A(I(1,i), I(2,i)) = a23(x1, h);
    A(I(2,i), I(1,i)) = a23(x1, h);
end
x1 = p(1, I(1,n));
x2 = p(1, I(2,n));
h = x2 - x1;
A(n, n) = ann(x1, h);
% b(n, 1) = bi(x1, x2, h);
b(n, 1) = bn(p(1, I(1,n)), h);
% 边值条件
x1 = p(1, I(1,1));
b(I(1,1), 1) = b(I(1,1), 1) - alpha.*a23(x1, h);
b(I(2,n-1), 1) = b(I(2,n-1), 1) + px(b_).*beta;
u = A\b;
c = cond(A);

dx = 1/1000;
x = a:dx:b_;
u1 = zeros(1, length(x));
for i=1:1:n-1
    for j=1:1:length(x)
        if p(1, I(1,i))<=x(1, j) && x(1, j)<=p(1, I(2,i)) % x_{i-1} <= x <= x_i
            h = p(1, I(2, i)) - p(1, I(1, i));
            u1(1, j) = u1(1, j) + ((x(1, j) - p(1, I(1, i))) / h) * u(I(1, i), 1);
        elseif p(1, I(1,i+1))<=x(1, j) && x(1, j)<=p(1, I(2,i+1)) % x_i <= x <= x_{i+1}
            h = p(I(2, i+1)) - p(I(1, i+1));
            u1(1, j) = u1(1, j) + (1 - (x(1, j) - p(I(1, i+1))) / h) * u(I(1, i), 1);
        end
    end
end
for j=1:1:length(x)
    if p(1, I(1,n))<=x(1, j) && x(1, j)<=p(1, I(2,n)) % x_{n-1} <= x <= x_n
        h = p(1, I(2, n)) - p(1, I(1, n));
        u1(1, j) = u1(1, j) + ((x(1, j) - p(1, I(1, n))) / h) * u(I(1, n), 1);
    end
end
for j=1:1:length(x)
    if p(1, I(1,1))<=x(1, j) && x(1, j)<=p(1, I(2,1)) % x_0 <= x <= x_1
        h = p(1, I(2, n)) - p(1, I(1, n));
        u1(1, j) = u1(1, j) + alpha * (1 - (x(1, j) - p(1, I(1, 1))) / h);
    end
end

u_ = sin(5.*x)+x.^3.*(2-x)+2;

% 绘制数值解和精确解的图像
figure(1);
plot(x, u1, 'b', 'LineWidth', 1.5);
hold on;
plot(x, u_, 'g', 'LineWidth', 1.5);
xlabel('x');
ylabel('u(x)');
legend('数值解', '精确解');
title('两点边值问题的有限元方法');
grid on;
hold off;
box off;

n=16时结果如下

在这里插入图片描述

n=32时结果如下

在这里插入图片描述

先求单元刚度矩阵

M a t l a b Matlab Matlab代码如下

% 先求单元刚度矩阵
n = ;
a = -1;
b_ = 2;

A = zeros(n, n);
b = zeros(n, 1);

p = zeros(1, n+1);
I = zeros(2, n);
for i=1:1:n+1
    p(1,i) = a + (b_-a)/n.*(i-1);
    if i<n+1
        I(1,i) = i;
        I(2,i) = i+1;
    end
end

u0 = @(x) sin(5*x)+x.^3.*(2-x)+2;
du = @(x) -3*(x - 2).*x.^2 - x.^3 + 5*cos(5*x);
alpha = u0(a);
beta = du(b_);
px = @(x) cos(x);
qx = @(x) x;
f = @(x) -((x - 2).*x.^3 - sin(5*x) - 2).*x + (6*(x - 2).*x + 6*x.^2 + 25*sin(5*x)).*cos(x) - (3*(x - 2).*x.^2 + x.^3 - 5*cos(5*x)).*sin(x);
a1 = @(x, h) integral(@(ksi) (px(x+h.*ksi)/h + h.*qx(x+h.*ksi).*(1-ksi).^2), 0, 1);
a23 = @(x, h) integral(@(ksi) (-px(x+h.*ksi)/h + h.*qx(x+h.*ksi).*ksi.*(1-ksi)),0, 1);
a4 = @(x, h) integral(@(ksi) (px(x+h.*ksi)/h + h.*qx(x+h.*ksi).*ksi.^2), 0, 1);
b1 = @(x, h) integral(@(ksi) (h.*f(x+h.*ksi).*(1-ksi)), 0, 1);
b2 = @(x, h) integral(@(ksi) (h.*f(x+h.*ksi).*ksi), 0, 1);

for i=1:1:n-1
    x1 = p(1, I(1,i+1));
    x2 = p(1, I(2,i+1));
    h = x2 - x1;
    a1_ = a1(x1, h);
    a2_ = a23(x1, h);
    a3_ = a23(x1, h);
    a4_ = a4(x1, h);
    A(I(1,i), I(1,i)) = A(I(1,i), I(1,i)) + a1_;
    A(I(1,i), I(2,i)) = A(I(1,i), I(2,i)) + a2_;
    A(I(2,i), I(1,i)) = A(I(2,i), I(1,i)) + a3_;
    A(I(2,i), I(2,i)) = A(I(2,i), I(2,i)) + a4_;
end
x1 = p(1, I(1,1));
x2 = p(1, I(2,1));
h = x2 - x1;
A(I(1,1), I(1,1)) = A(I(1,1), I(1,1)) + a4(x1, h);
for i=2:1:n
    x1 = p(1, I(1,i));
    x2 = p(1, I(2,i));
    h = x2 - x1;
    b1_ = b1(x1, h);
    b2_ = b2(x1, h);
    b(I(1,i-1), 1) = b(I(1,i-1), 1) + b1_;
    b(I(2,i-1), 1) = b(I(2,i-1), 1) + b2_;
end
x1 = p(1, I(1,1));
x2 = p(1, I(2,1));
h = x2 - x1;
b(I(1,1), 1) = b(I(1,1), 1) + b2(x1, h);
% 边值条件
b(I(1,1), 1) = b(I(1,1), 1) - alpha.*a23(x1, h);
b(I(2,n-1), 1) = b(I(2,n-1), 1) + px(b_).*beta;
u = A\b;
c = cond(A);

dx = 1/1000;
x = a:dx:b_;
u1 = zeros(1, length(x));
for i=1:1:n-1
    for j=1:1:length(x)
        if p(1, I(1,i))<=x(1, j) && x(1, j)<=p(1, I(2,i)) % x_{i-1} <= x <= x_i
            h = p(1, I(2, i)) - p(1, I(1, i));
            u1(1, j) = u1(1, j) + ((x(1, j) - p(1, I(1, i))) / h) * u(I(1, i), 1);
        elseif p(1, I(1,i+1))<=x(1, j) && x(1, j)<=p(1, I(2,i+1)) % x_i <= x <= x_{i+1}
            h = p(I(2, i+1)) - p(I(1, i+1));
            u1(1, j) = u1(1, j) + (1 - (x(1, j) - p(I(1, i+1))) / h) * u(I(1, i), 1);
        end
    end
end
for j=1:1:length(x)
    if p(1, I(1,n))<=x(1, j) && x(1, j)<=p(1, I(2,n)) % x_{n-1} <= x <= x_n
        h = p(1, I(2, n)) - p(1, I(1, n));
        u1(1, j) = u1(1, j) + ((x(1, j) - p(1, I(1, n))) / h) * u(I(1, n), 1);
    end
end
for j=1:1:length(x)
    if p(1, I(1,1))<=x(1, j) && x(1, j)<=p(1, I(2,1)) % x_0 <= x <= x_1
        h = p(1, I(2, n)) - p(1, I(1, n));
        u1(1, j) = u1(1, j) + alpha * (1 - (x(1, j) - p(1, I(1, 1))) / h);
    end
end

u_ = sin(5.*x)+x.^3.*(2-x)+2;

% 绘制数值解和精确解的图像
figure(1);
plot(x, u1, 'r', 'LineWidth', 1.5);
hold on;
plot(x, u_, 'b', 'LineWidth', 1.5);
xlabel('x');
ylabel('u(x)');
legend('数值解', '精确解');
title('两点边值问题的有限元方法');
grid on;
hold off;
box off;

n=16时结果如下

在这里插入图片描述

n=32时结果如下

在这里插入图片描述

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/684081.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

陶哲轩甩出调教GPT-4聊天记录,点击领取大佬的研究助理

量子位 | 公众号 QbitAI 鹅妹子嘤&#xff0c;天才数学家陶哲轩搞数学研究&#xff0c;已经离不开普通人手里的“数学菜鸡”GPT了&#xff01; 就在他最新解决的一个数学难题下面&#xff0c;陶哲轩明确指出自己“使用了GPT-4”&#xff0c;后者给他提出了一种可行的解决方法…

【FFmpeg实战】avformat_find_stream_info() 函数源码解析

转载自地址&#xff1a;https://cloud.tencent.com/developer/article/1873836 先来看一下 avformat_find_stream_info() 的头文件里的注释对该函数的介绍&#xff0c;本文我们基于 FFmpeg n4.2 版本的源码分析。 /*** Read packets of a media file to get stream informatio…

Apikit 自学日记:创建 API 文档

Apikit 中一共有5种创建API文档的方式&#xff1a; 新建API文档 导入API文档&#xff0c;详情可查看《导入、导出API文档》 从模板添加API文档&#xff0c;详情可查看《API文档模板》 自动生成API文档&#xff0c;详情可查看《自动生成API文档》 IDEA插件注释同步API文档 …

linux 在线安装 Redis

博主介绍&#xff1a; ✌博主从事应用安全和大数据领域&#xff0c;有8年研发经验&#xff0c;5年面试官经验&#xff0c;Java技术专家✌ Java知识图谱点击链接&#xff1a;体系化学习Java&#xff08;Java面试专题&#xff09; &#x1f495;&#x1f495; 感兴趣的同学可以收…

生成式AI掀起产业智能化新浪潮|爱分析报告

报告摘要 大模型支撑的生成式AI&#xff0c;让人类社会有望步入通用人工智能时代&#xff0c;拥有广阔的应用前景&#xff0c;有望赋能千行百业。当前生成式AI的落地整体处于初级阶段&#xff0c;不同模态的落地时间表差异明显&#xff0c;企业需求主要集中在数字化程度高、容…

地平线旭日x3派部署yolov8

地平线旭日x3派部署yolov8 总体流程1.导出onnx模型导出YOLOV8_onnxruntime.py验证onnxutils.py 2.在开发机转为bin模型2.1准备数据图片2.2转换必备的yaml文件2.3开始转换 3.开发机验证**quantized_model.onnx4.板子运行bin模型 资源链接 总体流程 1.导出onnx模型 导出 使用y…

03 | 事务隔离:为什么你改了我还看不见?

以下出自 《MySQL 实战 45 讲》 03 | 事务隔离&#xff1a;为什么你改了我还看不见&#xff1f; 隔离性与隔离级别 当数据库上有多个事务同时执行的时候&#xff0c;就可能出现脏读&#xff08;dirty read&#xff09;、不可重复读&#xff08;non-repeatable read&#xff0…

搜索功能全流程解析

在产品中一般会分布着大大小小的搜索&#xff0c;以便提升用户的信息获取效率和信息消费的能力。本文作者全流程角度&#xff0c;对搜索功能进行了讲解&#xff0c;并从搜索流程中寻找提升体验的触点&#xff0c;一起来看一下吧。 在产品中因多功能诉求和业务复杂性等因素&…

《Pytorch深度学习和图神经网络(卷 1)》学习笔记——第三章

学习基于如下书籍&#xff0c;仅供自己学习&#xff0c;用来记录回顾&#xff0c;非教程。 <PyTorch深度学习和图神经网络&#xff08;卷 1&#xff09;——基础知识>一书配套代码&#xff1a; https://github.com/aianaconda/pytorch-GNN-1st 百度网盘链接&#xff1a;…

vite优化

1.利用 rollup-plugin-analyzer 插件进行进行代码体积分析&#xff0c;从而优化你的代码。 根据项目体积分析&#xff0c;进行接下来的优化&#xff1a; &#xff08;一&#xff09;使用unplugin-vue-components插件按需加载antd vue 组件&#xff1a; 使用步骤 1、安装插件…

6.18 、Java初级:锁

1 同步锁 1.1 前言 经过前面多线程编程的学习,我们遇到了线程安全的相关问题,比如多线程售票情景下的超卖/重卖现象. 上节笔记点这里-进程与线程笔记 我们如何判断程序有没有可能出现线程安全问题,主要有以下三个条件: 在多线程程序中 有共享数据 多条语句操作共享数据 多…

GPT-4 的创造力全方位持平或碾压人类 | 一项最新研究发现

文章目录 一、前言二、主要内容三、总结 &#x1f349; CSDN 叶庭云&#xff1a;https://yetingyun.blog.csdn.net/ 一、前言 最近&#xff0c;一项有关 GPT-4 的创造力思维测试火了。来自蒙大拿大学和 UM Western 大学的研究团队发现&#xff0c;GPT-4 在 Torrance 创造性思维…

Sharding-JDBC之RangeShardingAlgorithm(范围分片算法)

目录 一、简介二、maven依赖三、数据库3.1、创建数据库3.2、创建表 四、配置&#xff08;二选一&#xff09;4.1、properties配置4.2、yml配置 五、范围分片算法六、实现6.1、实体层6.2、持久层6.3、服务层6.4、测试类6.4.1、保存订单数据6.4.2、根据时间范围查询订单 一、简介…

还在等待本地渲染?云渲染才是真的省时省心又省钱!

可能很多设计师会觉得本地渲染效果图或动画更灵活&#xff0c;而且没有什么额外的附加费用&#xff0c;但其实不然&#xff01;当你面对多个大型或紧急的项目时&#xff0c;本地渲染就“慌”了。 接下来我将全面对比“本地渲染”和“云渲染”&#xff0c;相信还在等待本地渲染…

黑客常用cmd命令(window版)

1、ping命令 ping命令是一个常用的网络工具&#xff0c;用来测试和诊断网络连接状况。通过发送ICMP&#xff08;Internet控制消息协议&#xff09;数据包到目标主机&#xff0c;并接收回复的数据包&#xff0c;可以测量目标主机的可达性、平均响应时间等指标。 在Windows操作…

前后端实现:行为验证码---文字点选

最近接到一个新的需求&#xff0c;由于客户是内网&#xff0c;你能使用腾讯的验证码了&#xff0c;需要改为前后端实现。 具体的代码已经提交git 项目效果图&#xff1a; 使用的技术栈&#xff1a;vitevue3ts git地址&#xff1a;https://github.com/susanliy/point_captcha…

TCP/IP协议是什么?

78. TCP/IP协议是什么&#xff1f; TCP/IP协议是一组用于互联网通信的网络协议&#xff0c;它定义了数据在网络中的传输方式和规则。作为前端工程师&#xff0c;了解TCP/IP协议对于理解网络通信原理和调试网络问题非常重要。本篇文章将介绍TCP/IP协议的概念、主要组成部分和工…

《程序喵》项目跨域问题解决思路

跨域问题&#xff1a;由于浏览器的 同源策略 限制&#xff0c;当一个请求url的协议、域名、端口号三者之间有任意一个与当前的url不同即为跨域。 同源策略是一种约定&#xff0c;它是浏览器中最核心也最基本的安全功能。同源策略会阻止一个域的 Javascript 脚本和另一个域的内…

举例说明梯度下降算法与最小二乘法的区别

梯度下降算法和最小二乘法都是用于求解线性回归问题中参数的优化方法。我们可以通过一个简单的例子来说明它们之间的区别。 假设我们有以下数据点&#xff1a;(1, 2)&#xff0c;(2, 3)&#xff0c;(3, 4)&#xff0c;(4, 5)&#xff0c;我们希望找到一条最佳拟合线 y wx b&a…

Android 中Looper机制详解

版本基于&#xff1a;Android R 0. 前言 在《Android 基于Handler 剖析消息机制》一文中&#xff0c;以 Handler 类为起点详细分析了异步通信&#xff0c;分析了Java 端 Handler 与Looper、MessageQueue、Message 之前的通信关系。 框架如下&#xff1a; 在Java 端的 Looper …