PDE с методом конечных разностей и итерацией Пикара в python - PullRequest
Новая
       12

PDE с методом конечных разностей и итерацией Пикара в python

0 голосов
/ 17 сентября 2018

Я хотел бы решить следующие PDE с питоном: 

\frac{\partial U}{\partial t}=\frac{\partial^2 U}{\partial x^2}-h(U,Q)\\
\frac{\partial Q}{\partial t}=h(U,Q)

, где h - заданная функция U и Q, со смешанными граничными условиями: 

Q(0,t)=1
\left[\frac{\partial U}{\partial x}\right]_{x=0}=-\alpha*(1-U(0,t))
Q(x\to\infty,t)=U(x\to\infty,t)=0

, где \ alphaэто заданная константа.Я использовал метод конечных разностей и итерацию Пикара.Все прошло гладко со следующим алгоритмом 

import numpy as np
import matplotlib.pyplot as plt
import os

## number of time steps 
n = 3200
## total time 
t = 1200.
## time increment 
dt =1.*t/n
## number of finite displacements 
k = 600
## maximal distance 
x = 200
## finite displacement 
dx = 1.*x/k

## initial conditions for U 
x0U = 1
u0 = 1.
## initial conditions for U have the form of a Heaviside step function 
u = np.zeros((n,k))
u[0,:x0U] = u0

## initial conditions for Q 
x0Q = 3
q0 = 1.
## the initial conditions for Q have the form of a Heaviside step function 
q = np.zeros((n,k))
q[0,:x0Q] = q0
## alpha
alpha = 10**(-1) 
## the following is the value for the boundary condition at infinity 
bond = 10.**(-25)

## define h
def h(y1,y2):
    return 1.*y1*(1.-y2)

## this is the matrix Mat used for Mat.u=b(u)
Mat = np.zeros((k,k))
for i in range(1,k-1):
    Mat[i,i+1]=1./(dx)**2
    Mat[i, i] = -1./dt - 2./(dx)**2
    Mat[i, i - 1] = 1./(dx)**2
Mat[0, 1] = 2./(dx)**2
Mat[0, 0] = -2./(dx)**2 - 2.*alpha/dx - 1./dt
Mat[k-1, k-1] = 1.
## define b(u)
b=np.zeros(k)
qQ = np.zeros(k)
uU = np.zeros(k)

## this is used to test whether the guess values are close enough to the results of the equation
epsilon = 10.**(-3.);
## the initial guess for u and q (called uG and qG) is the value in the previous time step
for i in range(1,n):
    qG = q[i-1]
    uG = u[i-1]
    if i%100==0:
        print   'uG ',uG[0]
        print   'qG ',qG[10]
## Picard's iteration starts here 
    check = 0;
    while check != 1:  
## find Q
        qQ = dt*h(uG, qG)+1.*q[i-1,:]
## update the constant term
        b = h(uG, qQ) -1.*u[i-1,:]/dt
        b[0] -= 2.*alpha/dx
        b[k-1] = bond
## solve the equation for U 
        uU = np.linalg.solve(Mat,b)
## check if the guess values are already close enough to the results of the equations
        if sum(abs(qQ - qG)) < epsilon and sum(abs(uU - uG)) < epsilon: 
            check = 1
## update Q and U
        else:   
            qG = qQ
            uG = uU
    if i%50==0: 
        print "time step ", i
## save the values for Q and U corresponding to time i
    q[i] = qQ
    u[i] = uU

fig1=plt.figure(1)
plt.plot(u[500])
plt.plot(u[1000])
plt.plot(u[1500])
plt.plot(u[2000]) 
fig5=plt.figure(2)
plt.plot(q[500])
plt.plot(q[1000])
plt.plot(q[1500])
plt.plot(q[2000]) 
plt.show()

Алгоритм вернет два графика, представляющих U и Q, как функцию от x для четырех различных значений t.Проблема возникает, когда я использую более сложную формулу для граничных условий, например 

\left[\frac{\partial U}{\partial x}\right]_{x=0}=-\frac{\alpha*(1-U(0,t))}{\alpha+(1-U(0,t))}

, где \ alpha - функция времени.В этих случаях бывают ситуации, когда U становится больше единицы, и это не должно происходить.Я хотел бы знать: 1) в целом, если мой алгоритм может быть улучшен и сделан быстрее 2), если это было хорошей идеей, использовать метод конечных разностей или есть более подходящие методы для этого конкретного случая (или, возможно, даже другие пакеты, которыесделай за меня работу) 3) если итерация Пикара была реализована правильно 4) почему я сталкиваюсь с этой численной проблемой для более сложных граничных условий и каковы общие приемы (если они есть), чтобы избежать проблемы

...