The focus of my research is on scientific computing and mathematical models' study. In particular, I study problems in tumor growth, blood coagulation, cell cycle control using both stochastic and deterministic approaches. The mathematical tools that I use include ODEs, PDEs, numerical algebraic geometry, uncertainty quantification, sensitivity analysis, bifurcation analysis, and computational methods.
Tumor growth is often modeled by PDEs with a free boundary where the changing shape of the tumor is of prime importance. Specifically, the models depend on values of the controlling parameter, called the tumor-aggressiveness factor, where bifurcations occur are nontrivial to compute. We have been using bifurcation analysis to prove that there exists a sequence of symmetric-breaking bifurcation branches of stationary solutions and using numerical algebraic geometry methods to compute the steady states of tumor growth models such as solid tumor, tumor with a necrotic core as well as fluid-like tumor.
Homotopy continuation to parameter \(\mu\)
In tumor models, it is interesting to find out whether it is possible
for the tumor to grow into other shapes besides spherical shape.
Homotopy continuation method combined with adaptive path tracking as well as multi-precision is used to compute nonradial solutions. Such procedure
automatically and systematically handles the singular system when the parameter approaches the bifurcation points. Numerical examples indicate that
the homotopy continuation method is efficient to manipulate parametric problems with bifurcation. Moreover, discretization of the free boundary by floating grids yields a polynomial system and reaches high order accuracy on the unconstructive grids.
W. Hao, J. D. Hauenstein, B. Hu Y. Liu, A. J. Sommese and Y.-T. Zhang, Continuation along bifurcation branches for a
tumor model with a necrotic core, Journal of Scientific Computing, to appear , (2012). [pdf]
W. Hao, J. D. Hauenstein, B. Hu T. McCoy and A. J. Sommese, Computing steady-state solutions for a free boundary problem modeling tumor growth by Stokes equation, Journal of Computational and Applied Mathematics, to appear, (2012). [pdf]
Blood coagulation model
Network of coagulation reactions
The structural integrity of a thrombus has important medical consequences, as fragments washed away from an unstable clot in a peripheral vein can embolize to the lungs with sometimes fatal results. The stability of the thrombus is influenced by the structural heterogeneity of the thrombus as the boundaries between discreet domains with different mechano-elastic properties are susceptible to fracture.
We develop algorithms based on numerical algebraic
geometry to compute steady state solutions for different initial
conditions. Because of rank-deficiency of the steady state system, a
numerical algebraic geometry algorithm is introduced to explore the
system and obtain the steady state solutions. Steady states which
can not be obtained by time marching method were discovered. We then
perform both variance decomposition based sensitivity analysis and
Morris design method to rank the significance of the 16 uncertain
reaction rates with respect to the total thrombin production.
Good
agreement on the sensitivity ranking of the reaction rates is
observed using both methods. Sparse grid probabilistic collocation
method (SGPCM) is employed to quantify how uncertainties of these
16 reaction rates influence total thrombin production. We also used the importance ranking of the corresponding
reaction rates obtained through sensitivity analysis study to
identify the most effective drug therapy and employ the model to
predict the
performance of the most effective drug therapy with confidence interval of
corresponding thrombin production using SGPCM.
The effect of platelet on blood coagulation is studied by developing
a distributed-Lagrange-multiplier/fictitious-domain(DLM/FD). Using
this method, the three-dimensional motion of a viscoelastic platelet
in a shear blood flow is simulated and compared with experiments on
tracking platelets in a blood chamber. The blood flow field is
modeled using an incompressible Navier-Stokes fluid solver; the
platelet is decoupled by Lagrange multiplier. We formulate potential
energy to model the forces of platelet, i.e., worm-like force,
bending force, surface force and volume force to allow accurate
representations of ''flipping" dynamics of platelets.
Related reference:
W. Hao, G. Lin, Z. Xu, A.J. Sommese and M. Alber, Analyzing stoichiometric regulation of blood coagulation using numerical algebraic geometry method and sensitivity analysis, Submitted, (2012). [pdf]
Numerical analysis/Numerical PDE
Bootstrapping method for computing multiple solutions of nonlinear PDEs
Discretization of systems of differential equations often leads to
systems of polynomials. Though using modern numerical codes for the
complete solution of polynomial systems can yield new solutions, the systems of polynomials (arising even from very sparse grids) are usually much too large for direct solution by
these codes. The realization underlying this article is that domain
decomposition gives excellent guidance on how to ``bootstrap'' from
the solutions of many small systems of polynomials to often large
numbers of solutions of a system of polynomials arising from a
discretization with a realistic grid. We introduces a new domain decomposition algorithm
which is called the bootstrapping method.
This method computes multiple solutions of nonlinear PDE systems
and is naturally parallelizable. Numerical examples for both one and two dimensional cases show the feasibility of algorithm.
Related references:
W. Hao, J. D. Hauenstein, B. Hu and A. J. Sommese, A domain decomposition algorithm for computing multiple steady states of differential equations, Submitted, (2011). [pdf]
Homotopy method based on WENO scheme for solving steady state problems of
hyperbolic conservation laws
Homotopy continuation is an
efficient tool originally designed for solving polynomial systems
via numerical algebraic geometry. Its efficiency relies on utilizing
adaptive stepsize and adaptive precision path tracking,
and endgames. In this paper, we apply homotopy
continuation to solve steady state problems of hyperbolic
conservation laws. The algorithm is based on a finite difference
discretization using the Lax-Friedrich numerical fluxes evaluated based
on the WENO scheme. Extensive numerical examples in both scalar and system
test problems in one and two dimensions demonstrate the efficiency.
Finite volume and finite difference methods are widely used in
hyperbolic conservation laws \[\mathbf{u}_t+\nabla\cdot
F(\mathbf{u})=g(\mathbf{u}).\] In both types of schemes, time
marching approach is expensive to compute steady states due to small
stepsize on time space. We are interested in developing homotopy method which is
based on WENO scheme and have a
faster convergence than the regular time marching approaches. Our
effort is restricted to steady state problems. We introduce a
homotopy function \[ H(u,\epsilon)=\big(\nabla\cdot
F(\mathbf{u})-g(\mathbf{u})-\epsilon \Delta
\mathbf{u}\big)(1-\epsilon)-\epsilon(\mathbf{u}-\mathbf{u}^0)\equiv0,\]
where \(\mathbf{u}^0\) is the initial condition, \(\epsilon\) is a
parameter between 0 and 1. In particular, the initial condition
automatically satisfies the homotopy function when \(\epsilon=1\), and
it becomes the steady state problem by setting
\(\epsilon=0\). Our homotopy method are tested on both one and two
dimensional examples and shows faster convergence than the time
marching approach. The figure below shows the result of
one dimensional Burger's equation.
Related references:
W. Hao, J. D. Hauenstein, A.J. Sommese, C.W. Shu, Z. Xu and Y. Zhang.Homotopy method for steady state problems on
hyperbolic conservation laws, Submitted, (2012) [pdf]
A package for computing multiplicity structures of nonlinear systems
Solving nonlinear systems of equations is one of the most basic problems in applied mathematics. It is well known that a solution of the system is multiple or
singular when the Jacobian is rank-deficient at the solution. In
other words, this particular solution is a converging point of a
few, a few dozens or even hundreds of solutions. Finding the
characteristics of the solution such as multiplicity is a natural
component of nonlinear system solving. However, the basic concept of
multiplicity of the solution and its identification are largely
unknown to numerical analysts and scientific computing
practitioners. None of the numerical analysis textbooks in all
levels we have seen elaborates the multiplicity beyond the single
univariate equations. We implemented a Matlab software, MULTIPLICITY, of a numerical algorithm for computing the multiplicity structure of a
nonlinear system at an isolated zero is presented. The software
incorporates a newly developed equation-by-equation strategy that
significantly improves the efficiency of the closedness subspace
algorithm and substantially reduces the storage requirement. As a
result, the algorithm and software can handle much larger nonlinear
systems and higher multiplicities than their predecessors.
Related references:
W. Hao, Andrew J. Sommese and Zhonggang Zeng, An algorithm and software for computing multiplicity structures at zeros of nonlinear systems, ACM Transactions on Mathematical Software, to appear, (2012). [pdf]
W. Hao, Andrew J. Sommese and Zhonggang Zeng, Matlab software for computing multiplicity structure of nonlinear systems, Available at here.
Parallel schemes for parabolic equations
Domain decomposition is a powerful tool for devising parallel
methods to solve time-dependent partial differential equations.
There is rich literature on domain decomposition finite difference
methods for solving parabolic equations on parallel computers. For
the non-overlapping domain decomposition methods, the explicit
nature of the calculation at the interface of sub-domain leads some
domain decomposition schemes to be conditionally stable, which
implies that they have to suffer from temporal step-size
restrictions. We developed parallel iterative method with unconditional stability and parallel finite difference scheme with high order accuracy as well as unconditional stability for solving parabolic equation. These schemes are based on domain decomposition method.
The numerical stability and convergence are derived in the
\(H^1\) norm in one dimensional case. Numerical results of two and three dimensions examine the stability, accuracy, and parallelism of the procedure.
Related references:
W. Hao and S. Zhu, A domain decomposition finite difference scheme with third-order accuracy and unconditional stability, submitted, (2012). [pdf]
I am also interested other interdisciplinary areas which numerical methods can be applied such as, numerical optimization on factorial design in Statistics; application of numerical algebraic geometry in the dorsal-ventral patterning of a zebrafish; mathematical modeling in gait control of the biped robot.
Gait control with data computed from a mathematical model we built up
Related references:
F. Sun, M.-Q. Liu and W. Hao, An algorithmic approach to finding factorial designs with generalized minimum aberration, Journal of Complexity, Volume 25, Issue 1, pp. 75-84, (2009). [pdf]
W. Hao and C. Liu, A high degree freedom mathematical model for the biped robot with the gait control, Mechanic in Engineering, Volume 28, No.6, pp. 69-72, (2006).
Discretization of systems of differential equations often leads to
systems of polynomials. Though using modern numerical codes for the
complete solution of polynomial systems can yield new solutions, the systems of polynomials (arising even from very sparse grids) are usually much too large for direct solution by
these codes. The realization underlying this article is that domain
decomposition gives excellent guidance on how to ``bootstrap'' from
the solutions of many small systems of polynomials to often large
numbers of solutions of a system of polynomials arising from a
discretization with a realistic grid. We introduces a new domain decomposition algorithm
which is called the bootstrapping method.
This method computes multiple solutions of nonlinear PDE systems
and is naturally parallelizable. Numerical examples for both one and two dimensional cases show the feasibility of algorithm.
