We investigate the time evolution and spatial variation of the stellar initial mass function (IMF) in star-forming disk galaxies by using chemodynamical simulations with an IMF model depending both on local densities and metallicities ([Fe/H]) of the interstellar medium (ISM). We find that the slope (α) of a power-law IMF (N(m)∝m -α) for stellar masses larger than 1 MOdot; evolves from the canonical Salpeter IMF (α ≈ 2.35) to be moderately top-heavy one (α ≈ 1.9) in the simulated disk galaxies with starbursts triggered by galaxy interaction. We also find that α in star-forming regions correlates with star formation rate densities (ΣSFR in units of MOdot; yr-1 kpc-2). Feedback effects of Type Ia and II supernovae are found to prevent IMFs from being too top-heavy (α <1.5). The simulation predicts α ≈ 0.23log ΣSFR + 1.7 for log ΣSFR ≥ -2 (i.e., more top-heavy in higher ΣSFR), which is reasonably consistent with corresponding recent observational results. The present study also predicts that inner regions of starburst disk galaxies have smaller α and thus are more top-heavy (dα/dR ∼ 0.07 kpc-1 for R ≤ 5 kpc). The predicted radial α gradient can be tested against future observational studies of the α variation in star-forming galaxies. © 2013. The American Astronomical Society. All rights reserved.