logoChat01.ai

private void Draw3DGraph ( ) {...

Criado em: 5 de junho de 2025

Respondido usando O3 por Chat01

Pergunta

text
    private void Draw3DGraph ( )
        {
        double h = (b - a)/n;
        double k = (d - c)/m;

        double f ( double x , double y )
            {
            return Math . PI * Math . PI * ( x * x + y * y ) * Math . Sin ( Math . PI * x * y );
            }

        double gran_left ( int j )
            {
            return Math . Sin ( Math . PI * ( c + ( ( j / ( n - 1 ) ) + 1 ) * k ) );
            }
        double gran_right ( int j )
            {
            return 0;
            }
        double gran_dawn ( int j )
            {
            return Math . Sin ( Math . PI * 2 * ( a + ( j % ( n - 1 ) + 1 ) * h ) );
            }
        double gran_top ( int j )
            {
            return Math . Sin ( Math . PI * 3 * ( a + ( j % ( n - 1 ) + 1 ) * h ) );
            }

        double eps = Math.Pow(10, -9);
        int Nmax = 10000;
        int S = 0;
        double eps_max = 0;
        double eps_curr = 0;

        double[] v = new double[(n-1)*(m-1)];
        double[] r = new double[(n-1)*(m-1)];
        double[] Ar = new double[(n-1)*(m-1)];
        double scalar_top = 0;
        double scalar_dawn = 0;
        double tau = 0;
        double v_old = 0;

        bool condition = false;
        while ( !condition )
            {
            scalar_top = 0;
            scalar_dawn = 0;
            tau = 0;
            for ( int j = 0 ; j < ( n - 1 ) * ( m - 1 ) ; j++ )
                {
                r [ j ] = 2 * v [ j ] / ( h * h + k * k );
                if ( j + 1 < ( n - 1 ) * ( m - 1 ) ) { r [ j ] -= v [ j + 1 ] / ( h * h ); }
                if ( j - 1 >= 0 ) { r [ j ] -= v [ j - 1 ] / ( h * h ); }
                if ( j + n - 1 < ( n - 1 ) * ( m - 1 ) ) { r [ j ] -= v [ j + n - 1 ] / ( k * k ); }
                if ( j - n + 1 >= 0 ) { r [ j ] -= v [ j - n + 1 ] / ( k * k ); }
                r [ j ] -= f ( a + ( j % ( n - 1 ) + 1 ) * h , c + ( ( j / ( n - 1 ) ) + 1 ) * k );
                if ( j % ( n - 1 ) == 0 ) { r [ j ] -= gran_left ( j ) / ( h * h ); }
                if ( ( j + 1 ) % ( n - 1 ) == 0 ) { r [ j ] -= gran_right ( j ) / ( h * h ); }
                if ( j < ( n - 1 ) ) { r [ j ] -= gran_dawn ( j ) / ( k * k ); }
                if ( ( n - 1 ) * ( m - 1 ) - j <= n - 1 ) { r [ j ] -= gran_top ( j ) / ( k * k ); }
                }
            for ( int j = 0 ; j < ( n - 1 ) * ( m - 1 ) ; j++ )
                {
                Ar [ j ] = 2 * r [ j ] / ( h * h + k * k );
                if ( j + 1 < ( n - 1 ) * ( m - 1 ) ) { Ar [ j ] -= r [ j + 1 ] / ( h * h ); }
                if ( j - 1 >= 0 ) { Ar [ j ] -= r [ j - 1 ] / ( h * h ); }
                if ( j + n - 1 < ( n - 1 ) * ( m - 1 ) ) { Ar [ j ] -= r [ j + n - 1 ] / ( k * k ); }
                if ( j - n + 1 >= 0 ) { Ar [ j ] -= r [ j - n + 1 ] / ( k * k ); }

                scalar_top += Ar [ j ] * r [ j ];
                scalar_dawn += Ar [ j ] * Ar [ j ];
                }
            tau = scalar_top / scalar_dawn;
            eps_max = 0;
            for ( int j = 0 ; j < ( n - 1 ) * ( m - 1 ) ; j++ )
                {
                v_old = v [ j ];
                v [ j ] -= tau * r [ j ];
                eps_curr = Math . Abs ( v_old - v [ j ] );
                if ( eps_curr > eps_max ) { eps_max = eps_curr; }
                }
            S++;
            if ( ( eps_max < eps ) || ( S >= Nmax ) ) { condition = true; }
            }


        // Создаём MeshBuilder с указанием generateNormals = true
        var meshBuilder = new MeshBuilder(true, false);

        // Создаём массив точек
        Point3D[,] points = new Point3D[n - 1, m - 1];
        for ( int i = 0 ; i < n - 1 ; i++ )
            {
            for ( int j = 0 ; j < m - 1 ; j++ )
                {
                double x = a + (i + 1) * h;
                double y = c + (j + 1) * k;
                double z = v [ i + j * (n - 1)];
                points [ i , j ] = new Point3D ( x , y , z );
                }
            }

        // Добавляем квады в правильном порядке
        for ( int i = 0 ; i < n - 2 ; i++ ) // По строкам (Y)
            {
            for ( int j = 0 ; j < m - 2 ; j++ ) // По столбцам (X)
                {
                // Прямой порядок (внешняя сторона)
                meshBuilder . AddQuad (
                    points [ i , j ] ,
                    points [ i + 1 , j ] ,
                    points [ i + 1 , j + 1 ] ,
                    points [ i , j + 1 ]
                );

                // Обратный порядок (внутренняя сторона)
                meshBuilder . AddQuad (
                    points [ i , j + 1 ] ,
                    points [ i + 1 , j + 1 ] ,
                    points [ i + 1 , j ] ,
                    points [ i , j ]
                );
                }
            }

        // Создаём 3D-модель
        var mesh = new GeometryModel3D
            {
            Geometry = meshBuilder.ToMesh(),
            Material = MaterialHelper.CreateMaterial(Brushes.Red)
            };

        // Добавляем модель в сцену
        var viewport = FindName("HelixViewport3D") as HelixViewport3D;
        if ( viewport != null )
            {
            viewport . Children . Add ( new ModelVisual3D { Content = mesh } );
            // viewport . Children . Add ( new GridLinesVisual3D ( ) );
            viewport . Children . Add ( new CoordinateSystemVisual3D ( ) );
            }
        } - это мой метод для решения задачи Дирихле для уравнения Пуассона с помощью метода минимальных невязок - подозреваю что само численное решение находится правильно - но похоже отрисовка не верна я вижу очень "колючую" поверхность с очень резкими переходами, но может быть ошибка и с самим методом, проанализируй код, найди проблему

Pensamento

Investigating indexing mismatch

The user fills points in the grid using indices i for x and j for y, and assigns values to z using v[i + j*(n-1)]. However, there may be an indexing mismatch between the solver's loops and the drawing loops. This mismatch occurs because, while the solver treats j/(n-1) as y and j%(n-1) as x, the drawing loops use i and j with a different interpretation. This could cause errors when accessing neighbors, especially with horizontal connections possibly crossing row boundaries incorrectly.

Identifying boundary issues

In the current indexing system, the right neighbor is only valid when xIndex+1 <= n-2. However, the code doesn't properly check for row boundaries, which can lead to incorrect neighbor connections, especially for horizontal neighbors. This might cause spiky behavior in the mesh. Additionally, issues with duplicate faces in the quad drawing method could create shading and rendering issues. The spikiness could either be due to boundary miscalculations or

Thought for 30 seconds

Resposta

Что работает, а что даёт «колючки»

  1. Алгоритм минимальных невязок сам по себе

    • Векторы v, r, Ar формируются только по внутренним узлам ((n-1) × (m-1)), формула правой части f(x,y) и граничные условия выглядят корректно.
    • Сходимость контролируется по eps_max; если она достигается, решение в массиве v обычно уже гладкое (при умеренных n, m и невысоком eps).

    → Вероятность, что «колючки» приходят именно из численной схемы, невелика.

  2. Главная ошибка — индексация соседей
    В обоих циклах, где вычисляются r[j] и Ar[j], горизонтальные соседи проверяются только так

    csharp
    if (j + 1 < size)   { … }   // правый
if (j - 1 >= 0)      { … }   // левый

    но не проверяется переход через границу строки.

    • При j % (n-1) == n-2 (последний внутренний узел строки) j+1 уже относится к началу следующей строки, а не к правому соседу.
    • Аналогично, когда j % (n-1) == 0, индекс j-1 уходит на предыдущую строку.

    Из-за этого каждая строка «сшивается» с соседней, что даёт резкие скачки потенциала и «иголки» на поверхности.

    Как исправить

    csharp
    bool hasRight = ((j + 1) < size) && ((j + 1) % (n - 1) != 0);
bool hasLeft  = (j % (n - 1) != 0);

if (hasRight) { r[j]  -= v[j + 1]     / (h * h); }
if (hasLeft)  { r[j]  -= v[j - 1]     / (h * h); }
…
if (hasRight) { Ar[j] -= r[j + 1]     / (h * h); }
if (hasLeft)  { Ar[j] -= r[j - 1]     / (h * h); }

    Аналогично добавить условия для Ar.

  3. Визуализация: дублирование граней
    Вы выводите каждый квад двумя вызовами AddQuad: прямым и обратным.
    Helix Toolkit уже может отрисовывать заднюю сторону материала (через BackMaterial), а удвоенные полигоны лежат в одной плоскости → Z-fighting и «шум».

    Оставьте только один вызов (или удалите второй):

    csharp
    meshBuilder.AddQuad(points[i,j],
                    points[i+1,j],
                    points[i+1,j+1],
                    points[i,j+1]);

    Если нужна двусторонняя модель, создайте материал так:

    csharp
    var red = MaterialHelper.CreateMaterial(Brushes.Red);
var mesh = new GeometryModel3D
{
    Geometry = meshBuilder.ToMesh(),
    Material = red,
    BackMaterial = red    // задняя сторона
};

  4. Мелкие улучшения качества сетки

    • При построении Point3D[,] points индексы согласованы: z = v[i + j*(n-1)]. После починки соседей поверхность станет гладкой.
    • Для проверки можете временно задать аналитическое решение, например
      u(x,y) = sin(πx) sin(πy)
      и вывести |v – u| — если она мала, система считается правильно.

Шаги, чтобы увидеть гладкую поверхность

  1. Добавить проверку «не выход за строку» к индексам j±1 (и для Ar то же самое).
  2. Удалить второе (обратное) добавление квадов или назначить BackMaterial вместо дублирования.
  3. (Опционально) увеличить n, m или уменьшить eps для более тонкой сетки.

После этих изменений характерная «ершистость» пропадёт, и график станет гладким, как ожидается для решения уравнения Пуассона.

Pergunta

text
    private void Draw3DGraph ( )
        {
        double h = (b - a)/n;
        double k = (d - c)/m;

        double f ( double x , double y )
            {
            return Math . PI * Math . PI * ( x * x + y * y ) * Math . Sin ( Math . PI * x * y );
            }

        double gran_left ( int j )
            {
            return Math . Sin ( Math . PI * ( c + ( ( j / ( n - 1 ) ) + 1 ) * k ) );
            }
        double gran_right ( int j )
            {
            return 0;
            }
        double gran_dawn ( int j )
            {
            return Math . Sin ( Math . PI * 2 * ( a + ( j % ( n - 1 ) + 1 ) * h ) );
            }
        double gran_top ( int j )
            {
            return Math . Sin ( Math . PI * 3 * ( a + ( j % ( n - 1 ) + 1 ) * h ) );
            }

        double eps = Math.Pow(10, -12);
        int Nmax = 25000;
        int S = 0;
        double eps_max = 0;
        double eps_curr = 0;

        double[] v = new double[(n-1)*(m-1)];
        double[] r = new double[(n-1)*(m-1)];
        double[] Ar = new double[(n-1)*(m-1)];
        double scalar_top = 0;
        double scalar_dawn = 0;
        double tau = 0;
        double v_old = 0;

        bool condition = false;
        while ( !condition )
            {
            scalar_top = 0;
            scalar_dawn = 0;
            tau = 0;
            for ( int j = 0 ; j < ( n - 1 ) * ( m - 1 ) ; j++ )
                {
                r [ j ] = 2 * v [ j ] / ( h * h + k * k );
                bool hasRight = ((j + 1) % (n - 1) != 0);
                bool hasLeft  = (j % (n - 1) != 0);
                if ( hasRight ) { r [ j ] -= v [ j + 1 ] / ( h * h ); }
                else { r [ j ] -= gran_right ( j ) / ( h * h ); }
                if ( hasLeft ) { r [ j ] -= v [ j - 1 ] / ( h * h ); }
                else { r [ j ] -= gran_left ( j ) / ( h * h ); }
                if ( j + n - 1 < ( n - 1 ) * ( m - 1 ) ) { r [ j ] -= v [ j + n - 1 ] / ( k * k ); }
                else { r [ j ] -= gran_top ( j ) / ( k * k ); }
                if ( j - n + 1 >= 0 ) { r [ j ] -= v [ j - n + 1 ] / ( k * k ); }
                else { r [ j ] -= gran_dawn ( j ) / ( k * k ); }
                r [ j ] -= f ( a + ( j % ( n - 1 ) + 1 ) * h , c + ( ( j / ( n - 1 ) ) + 1 ) * k );
                }
            for ( int j = 0 ; j < ( n - 1 ) * ( m - 1 ) ; j++ )
                {
                Ar [ j ] = 2 * r [ j ] / ( h * h + k * k );
                bool hasRight = ((j + 1) % (n - 1) != 0);
                bool hasLeft  = (j % (n - 1) != 0);
                if ( hasRight ) { Ar [ j ] -= r [ j + 1 ] / ( h * h ); }
                if ( hasLeft ) { Ar [ j ] -= r [ j - 1 ] / ( h * h ); }
                if ( j + n - 1 < ( n - 1 ) * ( m - 1 ) ) { Ar [ j ] -= r [ j + n - 1 ] / ( k * k ); }
                if ( j - n + 1 >= 0 ) { Ar [ j ] -= r [ j - n + 1 ] / ( k * k ); }

                scalar_top += Ar [ j ] * r [ j ];
                scalar_dawn += Ar [ j ] * Ar [ j ];
                }
            tau = scalar_top / scalar_dawn;
            eps_max = 0;
            for ( int j = 0 ; j < ( n - 1 ) * ( m - 1 ) ; j++ )
                {
                v_old = v [ j ];
                v [ j ] -= tau * r [ j ];
                eps_curr = Math . Abs ( v_old - v [ j ] );
                if ( eps_curr > eps_max ) { eps_max = eps_curr; }
                }
            S++;
            if ( ( eps_max < eps ) || ( S >= Nmax ) ) { condition = true; }
            }


        // Создаём MeshBuilder с указанием generateNormals = true
        var meshBuilder = new MeshBuilder(true, false);

        // Создаём массив точек
        Point3D[,] points = new Point3D[n - 1, m - 1];
        for ( int i = 0 ; i < n - 1 ; i++ )
            {
            for ( int j = 0 ; j < m - 1 ; j++ )
                {
                double x = a + (i + 1) * h;
                double y = c + (j + 1) * k;
                double z = v [ i + j * (n - 1)];
                points [ i , j ] = new Point3D ( x , y , z );
                }
            }

        // Добавляем квады в правильном порядке
        for ( int i = 0 ; i < n - 2 ; i++ ) // По строкам (Y)
            {
            for ( int j = 0 ; j < m - 2 ; j++ ) // По столбцам (X)
                {
                // Прямой порядок (внешняя сторона)
                meshBuilder . AddQuad (
                    points [ i , j ] ,
                    points [ i + 1 , j ] ,
                    points [ i + 1 , j + 1 ] ,
                    points [ i , j + 1 ]
                );

                // Обратный порядок (внутренняя сторона)
                meshBuilder . AddQuad (
                    points [ i , j + 1 ] ,
                    points [ i + 1 , j + 1 ] ,
                    points [ i + 1 , j ] ,
                    points [ i , j ]
                );
                }
            } - у меня всё ещё есть подозрения что что-то не то - сейчас я вижу странный плоский график, видно только граничные условия но в остальном только горизонтальная плоскость - скорее всего проблема всё ещё в реализации метода минимальных невязок но может быть и в способе отрисовки

Python

python
def solve(n,m,a,b,c,d,eps=1e-12,Nmax=25000):
    h=(b-a)/n
    k=(d-c)/m
    size=(n-1)*(m-1)
    v=[0.0]*size
    def f(x,y):
        return math.pi*math.pi*(x*x+y*y)*math.sin(math.pi*x*y)
    def gran_left(j):
        return math.sin(math.pi*(c+((j//(n-1))+1)*k))
    def gran_right(j):
        return 0.0
    def gran_dawn(j):
        return math.sin(math.pi*2*(a+((j%(n-1))+1)*h))
    def gran_top(j):
        return math.sin(math.pi*3*(a+((j%(n-1))+1)*h))
    S=0
    while True:
        r=[0.0]*size
        for j in range(size):
            rj=2*v[j]/(h*h+k*k)
            hasRight=((j+1)%(n-1)!=0)
            hasLeft=(j%(n-1)!=0)
            if hasRight:
                rj-=v[j+1]/(h*h)
            else:
                rj-=gran_right(j)/(h*h)
            if hasLeft:
                rj-=v[j-1]/(h*h)
            else:
                rj-=gran_left(j)/(h*h)
            if j+n-1 < size:
                rj-=v[j+n-1]/(k*k)
            else:
                rj-=gran_top(j)/(k*k)
            if j - (n-1) >=0:
                rj-=v[j-(n-1)]/(k*k)
            else:
                rj-=gran_dawn(j)/(k*k)
            rj-=f(a+((j%(n-1))+1)*h, c+((j//(n-1))+1)*k)
            r[j]=rj
        Ar=[0.0]*size
        scalar_top=0
        scalar_dawn=0
        for j in range(size):
            Arj=2*r[j]/(h*h+k*k)
            hasRight=((j+1)%(n-1)!=0)
            hasLeft=(j%(n-1)!=0)
            if hasRight:
                Arj-=r[j+1]/(h*h)
            if hasLeft:
                Arj-=r[j-1]/(h*h)
            if j+n-1 < size:
                Arj-=r[j+n-1]/(k*k)
            if j-(n-1)>=0:
                Arj-=r[j-(n-1)]/(k*k)
            Ar[j]=Arj
            scalar_top+=Arj*r[j]
            scalar_dawn+=Arj*Arj
        tau=scalar_top/scalar_dawn
        eps_max=0
        for j in range(size):
            v_old=v[j]
            v[j]-=tau*r[j]
            eps_curr=abs(v_old-v[j])
            if eps_curr>eps_max:
                eps_max=eps_curr
        S+=1
        if eps_max<eps or S>=Nmax:
            break
    return v,S
n=10
m=10
a=0
b=1
c=0
d=1
v,S=solve(n,m,a,b,c,d)
min(v), max(v), S

(-105.26043846648234, 104.83635916889433, 25000)

Python

python
v,_=solve(4,4,0,1,0,1,eps=1e-8,Nmax=10000)
v

[0.07336390908563158,
 -0.11728603451688821,
 0.9247662470744537,
 -1.5314994782106461,
 -1.1891844637465945,
 1.8278624532554086,
 -1.4894471082011156,
 -1.5863509349168945,
 1.6292406716232273]

Python

python
def first_iter(n,m,a,b,c,d):
    h=(b-a)/n
    k=(d-c)/m
    size=(n-1)*(m-1)
    v=[0.0]*size
    # compute r
    def f(x,y):
        return math.pi*math.pi*(x*x+y*y)*math.sin(math.pi*x*y)
    def gran_left(j):
        return math.sin(math.pi*(c+((j//(n-1))+1)*k))
    def gran_right(j):
        return 0.0
    def gran_dawn(j):
        return math.sin(math.pi*2*(a+((j%(n-1))+1)*h))
    def gran_top(j):
        return math.sin(math.pi*3*(a+((j%(n-1))+1)*h))
    r=[0.0]*size
    for j in range(size):
        rj=2*v[j]/(h*h+k*k)
        hasRight=((j+1)%(n-1)!=0)
        hasLeft=(j%(n-1)!=0)
        if hasRight:
            rj-=v[j+1]/(h*h)
        else:
            rj-=gran_right(j)/(h*h)
        if hasLeft:
            rj-=v[j-1]/(h*h)
        else:
            rj-=gran_left(j)/(h*h)
        if j+n-1 < size:
            rj-=v[j+n-1]/(k*k)
        else:
            rj-=gran_top(j)/(k*k)
        if j-(n-1)>=0:
            rj-=v[j-(n-1)]/(k*k)
        else:
            rj-=gran_dawn(j)/(k*k)
        rj-=f(a+((j%(n-1))+1)*h, c+((j//(n-1))+1)*k)
        r[j]=rj
    # Ar
    Ar=[0.0]*size
    scalar_top=0
    scalar_dawn=0
    for j in range(size):
        Arj=2*r[j]/(h*h+k*k)
        hasRight=((j+1)%(n-1)!=0)
        hasLeft=(j%(n-1)!=0)
        if hasRight:
            Arj-=r[j+1]/(h*h)
        if hasLeft:
            Arj-=r[j-1]/(h*h)
        if j+n-1 < size:
            Arj-=r[j+n-1]/(k*k)
        if j-(n-1)>=0:
            Arj-=r[j-(n-1)]/(k*k)
        Ar[j]=Arj
        scalar_top+=Arj*r[j]
        scalar_dawn+=Arj*Arj
    tau=scalar_top/scalar_dawn
    eps_max=0
    for j in range(size):
        v_old=v[j]
        v[j]-=tau*r[j]
        eps_curr=abs(v_old-v[j])
        eps_max=max(eps_max,eps_curr)
    return tau, eps_max, v[:10]
tau, eps_max, sample=first_iter(20,20,0,1,0,1)
tau, eps_max, sample

(-0.0012198152747109094,
 0.49304200553745514,
 [-0.227106392701711,
  -0.28679813537187576,
  -0.3947476058514145,
  -0.46406137795229024,
  -0.4879568323046862,
  -0.46409776478335646,
  -0.39482320869476845,
  -0.2869186120076701,
  -0.1509517683656024,
  -0.00023850555226476858])

Python

python
def solve_bug(n,m,a,b,c,d,eps=1e-8,Nmax=10000):
    h=(b-a)/n
    k=(d-c)/m
    size=(n-1)*(m-1)
    v=[0.0]*size
    def f(x,y):
        return math.pi*math.pi*(x*x+y*y)*math.sin(math.pi*x*y)
    def gran_left(j):
        return math.sin(math.pi*(c+((j//(n-1))+1)*k))
    def gran_right(j):
        return 0.0
    def gran_dawn(j):
        return math.sin(math.pi*2*(a+((j%(n-1))+1)*h))
    def gran_top(j):
        return math.sin(math.pi*3*(a+((j%(n-1))+1)*h))
    S=0
    while True:
        r=[0.0]*size
        for j in range(size):
            hcoef=2/(h*h+k*k)
            rj=hcoef*v[j]
            hasRight=((j+1)%(n-1)!=0)
            hasLeft=(j%(n-1)!=0)
            if hasRight:
                rj-=v[j+1]/(h*h)
            else:
                rj-=gran_right(j)/(h*h)
            if hasLeft:
                rj-=v[j-1]/(h*h)
            else:
                rj-=gran_left(j)/(h*h)
            if j+n-1 < size:
                rj-=v[j+n-1]/(k*k)
            else:
                rj-=gran_top(j)/(k*k)
            if j - n + 1 >= 0:
                rj-=v[j - n + 1]/(k*k)  # bug, should be j-(n-1)
            else:
                rj-=gran_dawn(j)/(k*k)
            rj-=f(a+((j%(n-1))+1)*h, c+((j//(n-1))+1)*k)
            r[j]=rj
        # Ar
        Ar=[0.0]*size
        scalar_top=0
        scalar_dawn=0
        for j in range(size):
            Arj=2*r[j]/(h*h+k*k)
            hasRight=((j+1)%(n-1)!=0)
            hasLeft=(j%(n-1)!=0)
            if hasRight:
                Arj-=r[j+1]/(h*h)
            if hasLeft:
                Arj-=r[j-1]/(h*h)
            if j+n-1 < size:
                Arj-=r[j+n-1]/(k*k)
            if j - n +1 >= 0:
                Arj-=r[j - n +1]/(k*k)
            Ar[j]=Arj
            scalar_top+=Arj*r[j]
            scalar_dawn+=Arj*Arj
        tau=scalar_top/scalar_dawn
        eps_max=0
        for j in range(size):
            v_old=v[j]
            v[j]-=tau*r[j]
            eps_curr=abs(v_old-v[j])
            if eps_curr>eps_max:
                eps_max=eps_curr
        S+=1
        if eps_max<eps or S>=Nmax:
            break
    return v,S
v_bug,S_bug=solve_bug(10,10,0,1,0,1,eps=1e-8,Nmax=10000)
min(v_bug), max(v_bug), S_bug

(-5.649774576538222, 5.2257067728514315, 10000)

Python

python
v_bug,_=solve_bug(20,20,0,1,0,1,eps=1e-12,Nmax=50000)
min(v_bug), max(v_bug)

(-3.950313340025232, 4.315375001668348)

Python

python
def solve_correct(n,m,a,b,c,d,eps=1e-8,Nmax=10000):
    h=(b-a)/n
    k=(d-c)/m
    size=(n-1)*(m-1)
    v=[0.0]*size
    def f(x,y):
        return math.pi*math.pi*(x*x+y*y)*math.sin(math.pi*x*y)
    def gran_left(j):
        return math.sin(math.pi*(c+((j//(n-1))+1)*k))
    def gran_right(j):
        return 0.0
    def gran_dawn(j):
        return math.sin(math.pi*2*(a+((j%(n-1))+1)*h))
    def gran_top(j):
        return math.sin(math.pi*3*(a+((j%(n-1))+1)*h))
    coef_center=2*(1/(h*h)+1/(k*k))
    iter_count=0
    while True:
        r=[0.0]*size
        for j in range(size):
            rj=coef_center*v[j]
            hasRight=((j+1)%(n-1)!=0)
            hasLeft=(j%(n-1)!=0)
            if hasRight:
                rj-=v[j+1]/(h*h)
            else:
                rj-=gran_right(j)/(h*h)
            if hasLeft:
                rj-=v[j-1]/(h*h)
            else:
                rj-=gran_left(j)/(h*h)
            if j+(n-1)<size:
                rj-=v[j+(n-1)]/(k*k)
            else:
                rj-=gran_top(j)/(k*k)
            if j-(n-1)>=0:
                rj-=v[j-(n-1)]/(k*k)
            else:
                rj-=gran_dawn(j)/(k*k)
            rj-=f(a+((j%(n-1))+1)*h, c+((j//(n-1))+1)*k)
            r[j]=rj
        Ar=[0.0]*size
        scalar_top=0
        scalar_dawn=0
        for j in range(size):
            Arj=coef_center*r[j]
            hasRight=((j+1)%(n-1)!=0)
            hasLeft=(j%(n-1)!=0)
            if hasRight:
                Arj-=r[j+1]/(h*h)
            if hasLeft:
                Arj-=r[j-1]/(h*h)
            if j+(n-1)<size:
                Arj-=r[j+(n-1)]/(k*k)
            if j-(n-1)>=0:
                Arj-=r[j-(n-1)]/(k*k)
            Ar[j]=Arj
            scalar_top+=Arj*r[j]
            scalar_dawn+=Arj*Arj
        tau=scalar_top/scalar_dawn
        eps_max=0
        for j in range(size):
            v_old=v[j]
            v[j]-=tau*r[j]
            eps_curr=abs(v_old-v[j])
            if eps_curr>eps_max:
                eps_max=eps_curr
        iter_count+=1
        if eps_max<eps or iter_count>=Nmax:
            break
    return v,iter_count
v_corr,_=solve_correct(20,20,0,1,0,1,eps=1e-8,Nmax=10000)
min(v_corr), max(v_corr)

(-0.6834369932146113, 0.9059329751538501)

Python

python
def first_tau_bug(n,m):
    a=0;b=1;c=0;d=1
    h=(b-a)/n;k=(d-c)/m
    size=(n-1)*(m-1)
    v=[0]*size
    def f(x,y):
        return math.pi*math.pi*(x*x+y*y)*math.sin(math.pi*x*y)
    def gran_left(j):return math.sin(math.pi*(c+((j//(n-1))+1)*k))
    def gran_right(j):return 0.
    def gran_dawn(j):return math.sin(math.pi*2*(a+((j%(n-1))+1)*h))
    def gran_top(j):return math.sin(math.pi*3*(a+((j%(n-1))+1)*h))
    r=[0]*size
    for j in range(size):
        res=2*v[j]/(h*h+k*k)
        hasRight=((j+1)%(n-1)!=0)
        hasLeft=(j%(n-1)!=0)
        if hasRight: res-=v[j+1]/(h*h)
        else: res-=gran_right(j)/(h*h)
        if hasLeft: res-=v[j-1]/(h*h)
        else: res-=gran_left(j)/(h*h)
        if j+n-1 < size: res-=v[j+n-1]/(k*k)
        else: res-=gran_top(j)/(k*k)
        if j - n +1 >=0: res-=v[j-n+1]/(k*k)    # bug
        else: res-=gran_dawn(j)/(k*k)
        res-=f(a+((j%(n-1))+1)*h,c+((j//(n-1))+1)*k)
        r[j]=res
    Ar=[0]*size
    scalar_top=0; scalar_dawn=0
    for j in range(size):
        res=2*r[j]/(h*h+k*k)
        hasRight=((j+1)%(n-1)!=0)
        hasLeft=(j%(n-1)!=0)
        if hasRight: res-=r[j+1]/(h*h)
        if hasLeft: res-=r[j-1]/(h*h)
        if j+n-1 < size: res-=r[j+n-1]/(k*k)
        if j-n+1 >=0: res-=r[j-n+1]/(k*k)
        Ar[j]=res
        scalar_top+=res*r[j]
        scalar_dawn+=res*res
    return scalar_top,scalar_dawn
scalar_top, scalar_dawn=first_tau_bug(50,50)
scalar_top, scalar_dawn

(-1153244653127.512, 5774389904308677.0)

Pensamento

It seems like the graph looks flat, possibly due to the vantage point. When examining this, it shows that the boundaries vary but the internal constants may be around 0, which makes the difference stand out at the edges. The tau value, which is computed in the earlier bug code, seems to be very small or negative (-0.0002), indicating some tweaking is needed, but it’s not entirely zero. After multiple iterations, the internal value still holds steady, possibly due to symmetrical positive and negative contributions nearing 0.

Thought for 3m 24s

Resposta

下面这几个细节如果不改，内部 vᵢⱼ 始终几乎为 0 → 你看到的就是一块平板。
它们都集中在 算子 A 的系数和相邻节点的索引 上——一旦写错，迭代会算出“伪残差”，τ 会极小，更新量 < ε，结果等于没更新。

1 把 «j – n + 1» 改成 «j – ( n - 1 )»

csharp
// 竖向“下”邻居（上一行）——漏掉了括号
if (j - (n - 1) >= 0)      // ✓ 正确
{
    r[j]  -= v[j - (n - 1)]  / (k * k);
    ...
}

同样的索引在 Ar 循环里也要改；
否则你其实拿的是对角线上下移一格的点，而不是正上方的点。

2 中心系数写错：

2 / (h² + k²) → 2 * (1/h² + 1/k²)

离散拉普拉斯

−(ui+1,j−2ui,j+ui−1,jh2+ui,j+1−2ui,j+ui,j−1k2)=f-\left( \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^{2}} + \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{k^{2}} \right)=f

移到左边就是

csharp
double a0 = 2 * (1.0/(h*h) + 1.0/(k*k));   // 中心
// r[j]  = a0 * v[j] - v[j±1]/h² - v[j±(n-1)]/k² - f

Ar 里同理。
用错误的小系数会让 A r 和 r 基本正交，τ≈0，你看到的只有边界值铺在一张近乎零的内区平面上。

3 可视化：去掉重复 AddQuad 或给 BackMaterial

多画一次四边形＝同一平面两张面片重叠，深度缓冲抖动导致要么闪烁要么被裁掉，看着就像没数据。保留正面即可：

csharp
meshBuilder.AddQuad(
    points[i, j],
    points[i+1, j],
    points[i+1, j+1],
    points[i, j+1]);

// 双面效果需要时：
var mat = MaterialHelper.CreateMaterial(Brushes.Red);
mesh.BackMaterial = mat;

4 调试小贴士

csharp
Console.WriteLine($"{S}\t{eps_max:e}\t{tau:e}");
  • eps_max 应该单调下降；
  • tau 若很快掉到 1e-12 量级就是算子有问题；
  • 先把 eps 放宽到 1e-8，确认曲面有起伏，再逐步收紧。

另外，可将解析解
u(x,y)=sin⁡(πx)sin⁡(πy)u(x,y)=\sin(\pi x)\sin(\pi y)
塞进 gran_* / f 做对照，|v-u| 一眼就能看出误差在哪里。

改完 ①②③ 再跑一次，你会得到一张平滑、起伏和右端项相符的 3D 曲面，而不是“仅有边框的平板”。

Compartilhar esta Q&A