偏微分方程数值解上机实验报告

matlab 偏微分 Freefem++ 计算

上机试验报告

第一题 求解 Helmholtz 方程的边值问题

题目

求解 Helmholtz 方程的边值问题：

$$\begin{eqnarray*}&&-\Delta u-k^2u=1,于 G=(0,1)\times(0,1),\\&&u=0, 于\quad\Gamma _1=\{x=0,0\leq y\leq 1\}\cup\{0\leq x\leq 1,y=1\}\\&&\frac{\partial u}{\partial n}=0, 于 \quad \Gamma _2=\{0\leq x\leq 1,y=0\}\cup\{x=1,0\leq y\leq1\} \end{eqnarray*}$$ 其中 $k=1,5,10,15,20$

分析

利用虚功原理求解 Helmholtz 方程的变分形式

$$\int\limits_G(-\Delta u)vdxdv+\int\limits_Gk^2uvdxdy=\int\limits_Gvdxdy$$ 因为 $$\Gamma _1=\{x=0,0\leq y\leq 1\}\cup\{0\leq x\leq 1,y=1\}，\\\Gamma _2=\{0\leq x\leq 1,y=0\}\cup\{x=1,0\leq y\leq1\}$$ 由格林第一公式 $$\begin{eqnarray*}\int\limits_G(-\Delta u)vdxdy&=&\int\limits_G\frac{\partial u}{\partial x}\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\frac{\partial v}{\partial y}dxdy-\int\limits_\Gamma\frac{\partial u}{\partial\textbf{n}}.vds\\&=&\int\limits_G\frac{\partial u}{\partial x}\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\frac{\partial v}{\partial y}dxdy\end{eqnarray*}$$

代码

使用 Freefem++ 编写程序实现使用有限元算法计算 Helmholtz 方程的边值问题,代码如下：

func z=0;
border a0(t=0,1){x=1;  y=t;}
border a1(t=0,1){x=1-t;y=1;}
border a2(t=0,1){x=0;  y=1-t;}
border a3(t=0,1){x=t;y=0;}
mesh Th = buildmesh(a0(150)+a1(150)+a2(150)+a3(150));
fespace Vh(Th,P2);
Vh phi,w,f=x*y;
solve Laplace(phi,w)=int2d(Th)(dx(phi)*dx(w) + dy(phi)*dy(w))
- int2d(Th)(f*w)+on(a1,phi=z)+on(a2,phi=z) ;
plot(phi,wait=true,ps="phi.eps");
plot(Th,wait=true, ps="membraneTh.eps");
savemesh(Th,"trace.msh");

结果

绘制的图形如下：

第二题 差分法求解 possion 问题

题目

用差分法求解如下 Laplace 方程 $(N=400)$

$$\left\{ \begin{array}{rcl} -\Delta u=xy,&&(x,y)\in \Omega=[0,1]\times [0,1] \\u=0,&&(x,y)\in\partial\Omega \end{array} \right.$$

分析

• 不妨记 $f(x,y) = xy, (x,y) \in \Omega = [0,1]\times [0,1]$ ,所以原方程可以写做 $-\Delta u=f(x,y)$ ，因为在边界出函数的值为零，所以是第一边值条件。
• 使用五点差分格式对方程进行差分
• 首先对函数的定义域进行剖分，根据 $N$ 的大小，可以确定网格的大小， 因此宽和高都为 $\frac{1}{N}$,以此也可以确定网点的位置。最终差分法会计算出函数在差分点的函数值，从而绘制出图形。
• 五点差分法的格式为 $$u_{ij}-\frac{1}{4}(u_{i-1,y}+u_{i,j-1}+u_{i+1,j}+u_{i,j+1})=\frac{h^2}{4}f_{ij}$$
• 利用 crs 压缩行存储以及五点差分法的格式可以获得差分方程系数矩阵的压缩存储格式以及值向量
• 利用基本的迭代法求解方程获得方程的解并绘制出图形
• 使用 freefem++ 通过变分法计算原方程的解并绘制出图形。与差分法绘制的图形进行对比。

代码

主函数 main

function [] = main()
N = input('lxease Input N:');
[value, cow, row, b]=getmatrix(N);
% for m=1:(N+1)^2
%     fprintf('\n')
%     for n = 1:(N+1)^2
%         %fprintf('%f',m)
%         temp = getvalue(m,n,value,cow,row);
%         fprintf('%f\t',temp)
%     end
% end
% b;
%% 使用基本迭代法求解方程组，（ E+ T）x = b.  then x = b-t*x
res  = solve_equation(value,cow,row,b);
%% now plot the shape
[y,x]=meshgrid(1:-1/N:0);
Z =flipud( reshape(res,N+1,N+1));
mesh(x,y,Z);
end

getmatrix 函数 主要用来获得系数矩阵的压缩行存储格式

function [value, cow, row,b] = getmatrix(N)
%% 使用addm 差值节点的地址信息。
[xsta, xend, ysta, yend] = deal(0, 1, 0, 1);%多变量同时赋值
xh = (xend-xsta)/N;
yh = (yend-ysta)/N;
addm = {};
for m = 1:N+1
    for n = 1:N+1
        xtemp = xsta + (m-1)*xh;
        ytemp = ysta + (n-1)*yh;
        addm(m,n) = {[xtemp, ytemp, xtemp*ytemp ],};
    end
end
%disp(addm);
%% 使用crs存储方法记录系数矩阵A
valuelen = (N*3-2)*N + N*(N-1)*2;
value = ones(valuelen,1);
cow = ones(N*N,1);
row = ones(valuelen,1);
tempadr = 1;%利用 tempadr 记录 value 数值变化的位置;
cowadr = 1;%利用 cowadr 记录 cow 列号变化的位置；
rowadr = 1 ;%利用 row 记录行变化的位置
b = zeros(N*N+2*N+1,1);
badr = 1;
%% 从第一行开始逐个记录系数矩阵每行的
for t = 1:N+1    
    if t==1
        for n = 1:N+1
            value(tempadr) = 1;
            tempadr = tempadr + 1;
            cow(cowadr) = tempadr;
            cowadr = cowadr +1;
            row(rowadr) = (t-1)*(N-1)+n;
            rowadr = rowadr + 1;
            b(badr) = 0;
            badr = badr + 1;
        end        
    else
        if t<=N               
            value(tempadr) = 1;
            tempadr = tempadr + 1;
            cow(cowadr) = tempadr;
            cowadr = cowadr +1;
            row(rowadr) = (t-1)*(N+1)+1;
            rowadr = rowadr + 1;
            b(badr) = 0;
            badr = badr + 1;
            for n=2:N
                value(tempadr:tempadr+4) = [-1/4  -1/4 1 -1/4 -1/4];
                tempadr = tempadr +5;
                cow(cowadr) = tempadr;
                cowadr = cowadr +1;
                row(rowadr:rowadr+4) = [(t-2)*(N+1)+n, (t-1)*(N+1)+n-1, (t-1)*(N+1)+n, (t-1)*(N+1)+n+1, t*(N+1)+n];
                rowadr = rowadr + 5;
                b(badr) = xh*xh*addm{t,n}(3);
                badr = badr + 1;
            end            
            value(tempadr) = 1;
            tempadr = tempadr + 1;
            cow(cowadr) = tempadr;
            cowadr = cowadr +1;
            row(rowadr) = t*(N+1);
            rowadr = rowadr + 1;
            b(badr) = 0;
            badr = badr + 1;                         
        else
            for n = 1:N+1
                value(tempadr) = 1;
                tempadr = tempadr + 1;
                cow(cowadr) = tempadr;
                cowadr = cowadr +1;
                row(rowadr) = (t-1)*(N+1)+n;
                rowadr = rowadr + 1;
                b(badr) = 0;
                badr = badr + 1;
            end             
        end
    end
end
end

solve_equation 主要使用基本迭代法求解方程组

function [ x ] = solve_equation( value, cow, row, b )
%This function is used to solve equation
%% get size of x;
lenth = length(b);
x0 = ones(lenth,1);
x = x0;
error = 1;
%% 方程的解为x
while error>=1e-5    
    for m=1:lenth
        temp = 0;
        if m == 1
            for t = 1:(cow(m)-1)
                if row(t) ~= m
                    n = row(t);
                    temp = temp + x0(n)*value(t);
                end
            end
        else
            for t = cow(m-1):(cow(m)-1)
                if row(t) ~= m
                    n = row(t);
                    temp = temp + x0(n)*value(t);
                end
            end
        end
        temp = b(m)- temp;
        x(m) = temp;       
    end  
   error = norm(x0-x,2);
   x0 = x;
end
end

getvalue 函数 检查系数矩阵是否正确

function[val] = getvalue(m,n,value,cow,row)
% 这个函数通过 crs 存储的 value , cow ,row 计算出矩阵第m行 n列的元素的值并返回
%%
lenth = length(cow);
if m>lenth || m<=0
    fprintf('Larger than the matrix')
else
    val = 0;
    if m==1
        for t=1:(cow(m)-1)
            if n == row(t)
                val =value(t);
                return;
            end
        end
    else
        val = 0;
        for t=cow(m-1):(cow(m)-1)
            if n == row(t)
                val =value(t);
                return;
            end
        end
    end
end

使用 freefem++ 通过变分法计算，代码如下：

real theta=4.*pi/3;
real a=2. ,b=1.;
func z=0;
border a0(t=0,1){x=1;  y=t;}
border a1(t=0,1){x=1-t;y=1;}
border a2(t=0,1){x=0;  y=1-t;}
border a3(t=0,1){x=t;y=0;}
mesh Th = buildmesh(a0(50)+a1(50)+a2(50)+a3(50));
fespace Vh(Th,P2);
Vh phi,w,f=x*y;
solve Laplace(phi,w)=int2d(Th)(dx(phi)*dx(w) + dy(phi)*dy(w))
- int2d(Th)(f*w)+on(a0,phi=z)+on(a1,phi=z)+on(a2,phi=z)+on(a3,phi=z) ;
plot(phi,wait=true);
plot(Th,wait=true, ps="membraneTh.eps");
savemesh(Th,"trace.msh");

结果

• 首先检查系数矩阵的正确性。输入的n 为2 矩阵的规模为 $9\times9$
• 程序运行结果所绘制的图形。
• 使用 freefem++ 绘制图形如下