2D基函数某个小分块矩阵的计算方法

托烦恼于北风

已于 2023-02-12 22:38:16 修改

阅读量82

点赞数

文章标签：矩阵线性代数

于 2023-02-12 20:31:23 首次发布

本文链接：https://blog.csdn.net/m0_63828036/article/details/128999037

版权

这段代码展示了如何利用基函数（如多项式）来构建4x4矩阵A，其中涉及到了二元函数的句柄定义和高斯积分的计算。代码中定义了四个基函数并计算了它们的不同组合，然后通过single2D_Gauss函数进行单元内部函数的积分。为了优化代码，作者提出了使用索引或元胞数组来减少重复性工作。

摘要由CSDN通过智能技术生成

基函数： $\phi {_{1}},,\phi {_{2}},\phi {_{3}},\phi {_{4}},\phi {_{5}},\phi {_{6}}$ .

$\phi _{1}=1, \phi _{2}=\frac{x-x_{i}}{\Delta x/2},\phi _{3}=\frac{y-y_{i}}{\Delta y/2},\\ ~~~~~\phi _{4}=\left (\frac{y-y_{i}}{\Delta y/2} \right )^2-\frac{1}{3},~~~~ \phi _{5}=\left (\frac{x-x_{i}}{\Delta x/2} \right )^2-\frac{1}{3}, \phi _{6}=\phi _{2}\phi _{3}.$

例如 $P^{1}$ 多项式，实际上需要算4x4矩阵A

$A=\begin{bmatrix} \int \phi _{_{1}}\phi _{_{1}}& \int \phi _{_{2}}\phi _{_{1}} & \int \phi _{_{3}}\phi _{_{1}} & \int \phi _{_{4}}\phi _{_{1}} \\ \int \phi _{_{1}}\phi _{_{2}} & \int \phi _{_{2}}\phi _{_{2}} & \int \phi _{_{3}}\phi _{_{2}} & \int \phi _{_{4}}\phi _{_{2}}\\ \int \phi _{_{1}}\phi _{_{3}} &\int \phi _{_{2}}\phi _{_{1}} & \int \phi _{_{3}}\phi _{_{3}} &\int \phi _{_{4}}\phi _{_{3}} \\ \int \phi _{_{1}}\phi _{_{4}} &\int \phi _{_{2}}\phi _{_{4}} &\int \phi _{_{3}}\phi _{_{4}} &\int \phi _{_{4}}\phi _{_{4}} \end{bmatrix}$ 。

下面是计算代码

a=0;
b=0.5;
c=0;
d=0.5;
dx=1/2;
dy=1/2;
xi=(a+b)/2;
yi=(c+d)/2;
A=zeros(4,4);
basis_1=@(x,y) 1;
basis_2=@(x,y) (x-xi)/(dx/2);
basis_3=@(x,y) (y-yi)/(dy/2);
basis_4=@(x,y) basis_2(x,y).* basis_3(x,y);



f11=@(x,y) basis_1(x,y).*basis_1(x,y);
f12=@(x,y) basis_2(x,y).*basis_1(x,y);
f13=@(x,y) basis_3(x,y).*basis_1(x,y);
f14=@(x,y) basis_4(x,y).*basis_1(x,y);
%2
f21=@(x,y) basis_1(x,y).*basis_2(x,y);
f22=@(x,y) basis_2(x,y).*basis_2(x,y);
f23=@(x,y) basis_3(x,y).*basis_2(x,y);
f24=@(x,y) basis_4(x,y).*basis_2(x,y);
%3
f31=@(x,y) basis_1(x,y).*basis_3(x,y);
f32=@(x,y) basis_2(x,y).*basis_3(x,y);
f33=@(x,y) basis_3(x,y).*basis_3(x,y);
f34=@(x,y) basis_4(x,y).*basis_3(x,y);
%4
f41=@(x,y) basis_1(x,y).*basis_4(x,y);
f42=@(x,y) basis_2(x,y).*basis_4(x,y);
f43=@(x,y) basis_3(x,y).*basis_4(x,y);
f44=@(x,y) basis_4(x,y).*basis_4(x,y);
%%%%%%%%%%%%%%%%%%%%%%%
A(1,1)=single2D_Gauss(a,b,c,d,f11,3);
A(1,2)=single2D_Gauss(a,b,c,d,f12,3);
A(1,3)=single2D_Gauss(a,b,c,d,f13,3);
A(1,4)=single2D_Gauss(a,b,c,d,f14,3);
%
A(2,1)=single2D_Gauss(a,b,c,d,f21,3);
A(2,2)=single2D_Gauss(a,b,c,d,f22,3);
A(2,3)=single2D_Gauss(a,b,c,d,f23,3);
A(2,4)=single2D_Gauss(a,b,c,d,f24,3);
%
A(3,1)=single2D_Gauss(a,b,c,d,f31,3);
A(3,2)=single2D_Gauss(a,b,c,d,f32,3);
A(3,3)=single2D_Gauss(a,b,c,d,f33,3);
A(3,4)=single2D_Gauss(a,b,c,d,f34,3);
%
A(4,1)=single2D_Gauss(a,b,c,d,f41,3);
A(4,2)=single2D_Gauss(a,b,c,d,f42,3);
A(4,3)=single2D_Gauss(a,b,c,d,f43,3);
A(4,4)= single2D_Gauss(a,b,c,d,f44,3);
A
%single2D_Gauss算的是单元内部函数f的积分（矩形单元的积分公式，建议自己编写！！）