How to create realistic shapes by interpolating two known shapes for facial blendshapes has not been investigated in the existing literature. In this paper, we propose a physics-based mathematical model and its analytical solutions to obtain more realistic facial shape changes. To this end, we first introduce the internal force of elastic beam bending into the equation of motion and integrate it with the constraints of two known shapes to develop the physics-based mathematical model represented with dynamic partial differential equations (PDEs). Second, we propose a unified mathematical expression of the external force represented with linear and various nonlinear time-dependent Fourier series, introduce it into the mathematical model to create linear and various nonlinear dynamic deformations of the curves defining a human face model, and derive analytical solutions of the mathematical model. Third, we evaluate the realism of the obtained analytical solutions in interpolating two known shapes to create new shape changes by comparing the shape changes calculated with the obtained analytical solutions and geometric linear interpolation to the ground-truth shape changes and conclude that among linear and various nonlinear PDE-based analytical solutions named as linear, quadratic, and cubic PDE-based interpolation, quadratic PDE-based interpolation creates the most realistic shape changes, which are more realistic than those obtained with the geometric linear interpolation. Finally, we use the quadratic PDE-based interpolation to develop a facial blendshape method and demonstrate that the proposed approach is more efficient than numerical physics-based facial blendshapes.