zwicker-group / py-pde

Python package for solving partial differential equations using finite differences.
https://py-pde.readthedocs.io
MIT License
402 stars 51 forks source link

why results is different from matlab #387

Closed lizhuang2511 closed 1 year ago

lizhuang2511 commented 1 year ago

Why use py-pde calculate results is different from use matlab Problem

Python code from math import pi

from pde import CartesianGrid, MemoryStorage, PDEBase, ScalarField, plot_kymograph,DiffusionPDE

class KortewegDeVriesPDE(PDEBase): """Korteweg-de Vries equation""" def init(self, k=54, diffusivity=[0.1,0.0025],bc="auto_periodic_neumann"): super().init() self.k = k self.ro=7272 self.c=420 self.a=self.k/(self.ro*self.c) self.q=7152 self.h1=39.6

self.diffusivity = diffusivity # spatial mobility

    bc_x_left = {"derivative":self.q/-self.k}
    bc_y = {"type": "mixed", "value": self.h1/self.k, "const":self.h1*22/self.k}
    self.bc = [bc_x_left,bc_y]  # boundary condition
    #self.bc = bc
def evolution_rate(self, state,t=0):
    """implement the python version of the evolution equation"""
    assert state.grid.dim == 1  # ensure the state is one-dimensional
    grad_x = state.laplace(self.bc)
    grad_x1=state.gradient(self.bc)[0]
    grad_x2=grad_x1.gradient(self.bc)[0]
    return grad_x*self.a

initialize the equation and the space

grid = CartesianGrid([[0, 10]],shape= [32], periodic=False)

grid = CartesianGrid([[0,0.055]],shape= [100])

state = ScalarField.from_expression(grid, "10")

state= ScalarField.random_normal(grid,20)

solve the equation and store the trajectory

storage = MemoryStorage() eq = KortewegDeVriesPDE() eq.solve(state, t_range=[0,100],tracker=storage.tracker(0.5))

plot the trajectory as a space-time plot

plot_kymograph(storage) python results

matlab code k=52 p=7272 c=420 a=k/p/c dd=1/a q=7152 ccc=q/k h=39.6 x = linspace(0,0.055,23); t = linspace(0,300,100); m = 0; sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t); u = sol(:,:,1); surf(x,t,u) title('Numerical solution computed with 20 mesh points') xlabel('Distance x') ylabel('Time t') function [c,f,s] = pdex1pde(x,t,u,dudx) c = 5.873538461538461e+04; f = dudx; s = 0 end

function u0 = pdex1ic(x) u0 = 20; end

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) pl = 1.375384615384616e+02; ql = 1; pr = -39.6xr+39.620; qr = -52; end matlab results

lizhuang2511 commented 1 year ago

![Uploading Screenshot_20230405_214556_com.foxit.mobile.pdf.lite.jpg…]()

david-zwicker commented 1 year ago

It's unfortunately not clear to me what you're asking or reporting. Please improve the report with a clear description of the problem, the expected behavior, and preferably code that can be run independently.