快速扫描算法求解Eikonal方程咨询

Question

我在使用快速扫描法（Fast sweeping method）求解程函方程时，发现写出的程序虽然能够收敛，但是数值解与真解之间的误差很大，请各位帮忙看看，以下是我的代码。其中使用的算例的是

参考的文章为：Zhao, H, K., 2005. A Fast sweeping method for Eikonal equations

%fast sweeping for Eikonal
function [phi_new, phi_old] = Eikonal(Nx, Ny, r, x, y, dx, dy, center, ghost_point, phi,c)
%c=1
d = @(x, y) sqrt((x - center(1))^2 + (y - center(2))^2); % 距离函数
h = dx;
phi_old = phi;
phi_new = phi;
% 快速扫掠法求解Eikonal方程 |∇phi_0| = c(x,y,T_end)
for i = 1 : 1 : Nx
    for j = 1 :1 : Ny         
        if i==1
            min_x = min(ghost_point, phi(i+1,j));
        elseif i==Nx
            min_x = min(phi(i-1,j),ghost_point);
        else
            min_x = min(phi(i-1,j),phi(i+1,j));
        end
        if j==1
            min_y = min(ghost_point, phi(i,j+1));
        elseif j==Ny
            min_y = min(phi(i,j-1),ghost_point);
        else
            min_y = min(phi(i,j+1), phi(i,j-1));
        end
        if abs(min_x-min_y)>=(c(i,j))*h
            tmp = min(min_x,min_y) + (c(i,j)) * h;
        else
            tmp = 0.5*(min_x+min_y+sqrt(abs(2*(c(i,j)^2)*(h*h)-(min_x-min_y)*(min_x-min_y))));
        end
        phi(i,j) = min(tmp,phi(i,j));
        phi_new(i,j) = phi(i,j);
        if x(i)==0 && y(j) ==0
            phi(i,j) = phi_old(i,j);
            phi_new(i,j) = phi(i,j);
        end
    end
end
for i = Nx : -1 : 1
    for j = 1 :1 : Ny
        if i==1
            min_x = min(ghost_point, phi(i+1,j));
        elseif i==Nx
            min_x = min(phi(i-1,j),ghost_point);
        else
            min_x = min(phi(i-1,j),phi(i+1,j));
        end
        if j==1
            min_y = min(ghost_point, phi(i,j+1));
        elseif j==Ny
            min_y = min(phi(i,j-1),ghost_point);
        else
            min_y = min(phi(i,j+1), phi(i,j-1));
        end
        if abs(min_x-min_y)>=(c(i,j))*h
            tmp = min(min_x,min_y) + (c(i,j)) * h;
        else
            tmp = 0.5*(min_x+min_y+sqrt(abs(2*(c(i,j)^2)*(h*h)-(min_x-min_y)*(min_x-min_y))));
        end
        phi(i,j) = min(phi(i,j),tmp);
        if x(i)==0 && y(j) ==0
            phi(i,j) = phi_old(i,j);
        end
    end
end
for i = Nx : -1 : 1
    for j = Ny :-1 : 1
        if i==1
            min_x = min(ghost_point, phi(i+1,j));
        elseif i==Nx
            min_x = min(phi(i-1,j),ghost_point);
        else
            min_x = min(phi(i-1,j),phi(i+1,j));
        end
        if j==1
            min_y = min(ghost_point, phi(i,j+1));
        elseif j==Ny
            min_y = min(phi(i,j-1),ghost_point);
        else
            min_y = min(phi(i,j+1), phi(i,j-1));
        end
        if abs(min_x-min_y)>=(c(i,j))*h
            tmp = min(min_x,min_y) + (c(i,j)) * h;
        else
            tmp = 0.5*(min_x+min_y+sqrt(abs(2*(c(i,j)^2)*(h*h)-(min_x-min_y)*(min_x-min_y))));
        end
        phi(i,j) = min(tmp,phi(i,j));
        if x(i)==0 && y(j) ==0
            phi(i,j) = phi_old(i,j);
        end
    end
end
for i = 1 : 1 : Nx
    for j = Ny :-1 : 1
        if i==1
            min_x = min(ghost_point, phi(i+1,j));
        elseif i==Nx
            min_x = min(phi(i-1,j),ghost_point);
        else
            min_x = min(phi(i-1,j),phi(i+1,j));
        end
        if j==1
            min_y = min(ghost_point, phi(i,j+1));
        elseif j==Ny
            min_y = min(phi(i,j-1),ghost_point);
        else
            min_y = min(phi(i,j+1), phi(i,j-1));
        end
        if abs(min_x-min_y)>=(c(i,j))*h
            tmp = min(min_x,min_y) + (c(i,j)) * h;
        else
            tmp = 0.5*(min_x+min_y+sqrt(abs(2*(c(i,j)^2)*(h*h)-(min_x-min_y)*(min_x-min_y))));
        end
        phi(i,j) = min(tmp,phi(i,j));
        if x(i)==0 && y(j) ==0
            phi(i,j) = phi_old(i,j);
        end
    end
end
end

clc; clear; close all;
h=[0.02,0.01,0.005,0.0025,0.002];
x = (-0.5:h(3):0.5);
y = (-0.5:h(3):0.5);
Nx = length(x);
Ny = length(y);
[X, Y] = meshgrid(x, y);
X = X';
Y = Y';
center = [0,0];
d = @(x, y) sqrt((x - center(1))^2 + (y - center(2))^2);
c = @(x,y) 1+exp(-40*((x+0.2)^2+(y+0.2)^2))-exp(-40*((x-0.2)^2+(y-0.2)^2)); %程函方程中的f
true_phi = @(x,y) abs(exp((sqrt(x^2+y^2)-0.2)*(0.01*x+x*y+y^2+0.4))-1);
ghost_point = 1e12;
phi = ones(Nx,Ny)*1e6;
c_used = zeros(Nx,Ny);
true_values = zeros(Nx,Ny);
r=0.015;
%初始化矩阵
for i = 1 : Nx
    for j =1 : Ny
        if x(i)==0 && y(j) ==0
            phi(i,j) = 0;
        end
        c_used(i,j) = c(x(i),y(j));
        true_values(i,j) = true_phi(x(i),y(j));
    end
end
%c_used = c(X,Y);
delta = 1;
k=1;
while delta>1e-9
    [phi_new,phi_old] = Eikonal(Nx,Ny, r,x,y,h(3),h(3),center,ghost_point,phi,c_used);
    delta = max(max(phi_new - phi_old));
    phi = phi_new;
    k = k+1;
end
error = max(max(phi-true_values));

以下是我的求解结果

以下为原文中的算例结果

Answer 1

以下参考免费微信小程序：皆我百晓生 作答，并由阿齐Archie整理回复，若还有疑问可评论或留言。

在您提供的代码中，您正在使用快速扫描法（Fast Sweeping Method）来求解Eikonal方程。Eikonal方程是一个非线性偏微分方程，通常用于计算最短路径或波前传播。在您的代码中，您尝试求解的是一个具有特定速度场的Eikonal方程。

根据您提供的代码和问题描述，以下是一些可能导致数值解与真实解之间误差较大的因素，以及相应的解决方案：

1. 初始化条件：

   • 确保初始条件phi在源点（即u(0,0)=0）附近设置为0，其他地方设置为一个足够大的值。这样可以确保波前从源点开始传播。

2. 边界条件：

   • 检查边界条件是否正确实现。在您的代码中，您使用了ghost_point来模拟边界条件，确保这是正确处理的。

3. 数值稳定性：

   • 检查时间步长（在这种情况下是空间分辨率dx和dy）是否足够小，以确保数值稳定性。如果步长太大，可能会导致数值解的不稳定。

4. 速度场c(x,y)：

   • 确保速度场c(x,y)在所有点上都是正的，因为Eikonal方程要求速度不能为零或负。

5. 迭代收敛：

   • 确保迭代足够多次，直到解收敛。您的代码中使用了delta来检查收敛性，但可能需要调整阈值或增加迭代次数。

6. 精度问题：

   • MATLAB中的数值计算可能会受到浮点精度的限制。检查是否有必要使用更高精度的数据类型，如vpa（符号计算）或double（双精度浮点数）。

7. 代码逻辑：

   • 仔细检查代码逻辑，确保在更新phi值时没有逻辑错误。

8. 真实解的验证：

   • 确保您的真实解（true_phi）是正确的，并且与您的数值方法和离散化方案相匹配。

针对您的代码，以下是一些具体的建议：

• 在更新phi值时，您使用了min函数来选择最小的值。这可能不是最佳选择，因为它可能导致解的下界被低估。相反，您应该使用正确的公式来更新phi值，该公式考虑了速度场c(x,y)和距离函数d(x,y)。

• 您可能需要重新审视您的c(x,y)函数，确保它在整个域内都是正的，并且符合物理意义。

• 您的true_phi函数看起来非常复杂，确保它正确地代表了真实的解析解。

• 您的迭代停止条件是基于delta的，但您可能需要更仔细地检查这个值是否真的表示了解的收敛。

最后，您可能需要对您的代码进行调试，以确保每个步骤都按照预期工作。您可以通过绘制中间步骤的phi值来可视化波前的传播，并检查它是否符合预期的行为。如果可能，您还可以与其他数值方法（如有限差分法或有限元法）的结果进行比较，以验证您的解。