This is a problem on flow around cylinder using total potential formulation with its analytical solution in Python © 2025 Shibo Peng

import numpy as np
import pandas as pd
import matpl.otlib.pyplot as plt
from mpl_toolkits import mplot3d
import warnings
warnings.filterwarnings('ignore')

def dip2p(x1,y1,x2,y2,xp,yp): #find subtended angle
    vec1 = [xp-x1,yp-y1]
    vec2 = [xp-x2,yp-y2]
    alpha = np.dot(vec1,vec2)/(np.sqrt(np.dot(vec1,vec1))*np.sqrt(np.dot(vec2,vec2)))
    return np.arccos(alpha)/2/np.pi

def poly_ellipse(A,B,n): #generate ellipse or cicle node points
    dt = 2*np.pi/n
    polyx = np.zeros(n)
    polyy = np.zeros(n)
    count = 0
    for i in np.arange(0,n):
        polyx[i] = A*np.cos(count*dt)
        polyy[i] = B*np.sin(count*dt)
        count += 1
        
    return polyx, polyy
def buildmatrix(xs,xc,ys,yc): #Aij builder
    n = np.size(xc)
    A = np.zeros((n,n))
    for k in np.arange(0,n):
            for i in np.arange(0,n):
                if i != k:
                    if i != n-1:
                        A[k,i] = dip2p(xs[i],ys[i],xs[i+1],ys[i+1],xc[k],yc[k])
                    else:
                        A[k,i] = dip2p(xs[i],ys[i],xs[0],ys[0],xc[k],yc[k])
                else: 
                    A[k,i] = 0.5
    return A
def control_points(xs,ys): #find control points
    n = xs.size
    xc = np.zeros(n)
    yc = np.zeros(n)
    for i in np.arange(0,n):
        if i != n-1:
            xc[i] = (xs[i]+xs[i+1])/2
            yc[i] = (ys[i]+ys[i+1])/2
        else: 
            xc[i] = (xs[i]+xs[0])/2
            yc[i] = (ys[i]+ys[0])/2
    return xc,yc

def solver(input,n,BA): #matrix solver
    dt = 2*np.pi/n
    xs = poly_ellipse(1,BA,n)[0]
    ys = poly_ellipse(1,BA,n)[1]
    xc = control_points(xs,ys)[0]
    yc = control_points(xs,ys)[1]
    A = buildmatrix(xs,xc,ys,yc)
    if input == 0:
        Phi_a = xc + xc/(xc**2+yc**2)
    else:
        Phi_a = (1+BA)*xc
    
    Phi_i = xc[np.newaxis].T
    Phi = np.dot(np.linalg.inv(A),Phi_i)

    angle = np.zeros(yc.size)
    for i in np.arange(0,yc.size):
        angle[i] = 360/2/np.pi*np.arctan2(yc[i],xc[i])
        if angle[i] < 0:
            angle[i] = angle[i] + 360
    return Phi_a, Phi, angle, xc, yc

def error(Phi_a,Phi): #find the total error 
    n = Phi_a.size
    err = 0
    for i in np.arange(0,n):
        if err <= abs(Phi_a[i]-Phi[i]):
            err = abs(Phi_a[i]-Phi[i])
            
    return err

def error_q(Phi_a,Phi): #find the error at pi/4
    n = Phi_a.size
    err = abs(Phi_a[int((n/4+1)/2)]-Phi[int((n/4+1)/2)])/abs(Phi_a[int((n/4+1)/2)])
    return err