first task of basic mathematical modeling

characteristic.py

# -*- coding: utf-8 -*-


import numpy as np
import matplotlib.pyplot as plt


def CharacteristicX(x0, t):
    a = 2+np.cos(np.cos(np.pi*x0/2))
    b = 1+(1+2*np.cos(np.pi*x0/2)+np.sin(np.cos(np.pi*x0/2)))**2
    return x0-a*t/b


def CharacteristicT(t0, t):
    a = 2 + np.cos(1+0.5*np.arctan(t0))
    b = 1 + (3+np.arctan(t0)+np.sin(1+0.5*np.arctan(t0)))**2
    return -a*(t-t0)/b
    

def BuildGraph(dataX, dataY, init, flg):
    s = "x0 = "
    if flg:
        s = "t0 = "
    graph = plt.figure()
    ax = graph.add_subplot(111)
    i = 0
    for m in dataY:
        ax.plot(dataX, m, label = s + str(init[i]))
        i += 1
    ax.set_ylim(-1, 0)
    ax.set_xlim(0, 6)
    plt.legend()
    ax.set_xlabel('t')
    ax.set_ylabel('x')
    plt.show()


def main():
    #t0 = 0
    mas = []
    T = np.arange(0, 7, 0.5)
    X = [-1, -0.8, -0.6, -0.4, -0.2, 0]
    for x0 in X:
        temp = []
        for t in T:
            temp.append(CharacteristicX(x0, t))
        mas.append(temp)
    BuildGraph(T, mas, X, 0)
    
    #x0 = 0
    mas = []
    T = np.arange(0, 7, 0.5)
    T0 = range(0, 7, 1)
    for t0 in T0:
        temp = []
        for t in T:
            temp.append(CharacteristicT(t0, t))
        mas.append(temp)
    BuildGraph(T, mas, T0, 1)
    
    
if __name__ == "__main__" :
    main()

main.py

# -*- coding: utf-8 -*-


import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


class Grid:

    def __init__(self, x0, x1, t0, t1, Nx, Nt):
        self.x0, self.x1 = x0, x1
        self.t0, self.t1 = t0, t1
        self.Nx, self.Nt = Nx, Nt
        self.h = (x1-x0)/Nx
        self.tau = (t1-t0)/Nt
        self.eps = 10**-6
        
        self.grid = np.zeros((Nx, Nt), dtype=float)
        self.dataX = np.linspace(x0, x1, Nx)
        self.dataT = np.linspace(t0, t1, Nt)
    
    
    def InitDataX(self, x):
        return np.cos(np.pi*x/2)
        
        
    def InitDataT(self, t):
        return 1+0.5*np.arctan(self.dataT)
        
        
    def InitGrid(self):
        self.grid[:, 0] = Grid.InitDataX(self, self.dataX)
        self.grid[0, :] = Grid.InitDataT(self, self.dataT)

    
    def F(self, x):
        return -np.arctan(1+2*x+np.sin(x))
    
    
    def FGrid(self, u11, u12, u21, u22):
        a = (u12-u11-u21+u22)/(2*self.tau)
        b = (Grid.F(self, u21)-Grid.F(self, u11)-Grid.F(self, u12)+Grid.F(self, u22)) / (2*self.h)
        return a+b
    
    
    def CoefficientC(self, x):
        return -(2+np.cos(x))/(1+(1+2*x+np.sin(x))**2)
        
    
    def DFGrid(self, x):
        return 1/(2*self.tau) + Grid.CoefficientC(self, x)/(2*self.h)
    
    
    def Newton(self, x_i, t_i):
        delta = self.eps+1
        u22_0 = self.grid[x_i+1, t_i+1]
        u11 = self.grid[x_i, t_i]
        u12 = self.grid[x_i, t_i+1]
        u21 = self.grid[x_i+1, t_i]

        while delta > self.eps:
            u22_1 = u22_0 - Grid.FGrid(self, u11, u12, u21, u22_0)/Grid.DFGrid(self, u22_0)
            delta = abs(u22_0-u22_1)
            u22_0 = u22_1
        
        return u22_0
        
    
    def FillGrid(self):
        for t_i in range(self.Nt-1):
            for x_i in range(self.Nx-1):
                self.grid[x_i+1, t_i+1] = Grid.Newton(self, x_i, t_i)
    
    
    def BuildGrpah(self):
        fig = plt.figure()
        ax = fig.add_subplot(111, projection='3d')
        self.dataX = self.dataX[::-1]
        self.dataX, self.dataT = np.meshgrid(self.dataX, self.dataT)
        surf = ax.plot_surface(self.dataT, self.dataX, self.grid, cmap='magma')
        ax.set_xlabel('t')
        ax.set_ylabel('x')
        ax.set_zlabel('u')
        plt.show()
        
        
    def FullSolution(self):
        Grid.InitGrid(self)
        Grid.FillGrid(self)
        Grid.BuildGrpah(self)

        
def main():
    solution = Grid(0, -1, 0, 6, 200, 200)
    solution.FullSolution()

    
if __name__ == "__main__" :
    main()