I am trying to convert the following Python example to Julia:
""" Tutorial example showing how to use the implicit solver RADAU. It simulates a falling mass. """ import numpy as np import pylab as plt from assimulo.problem import Implicit_Problem # Imports the problem formulation from Assimulo from assimulo.solvers import Radau5DAE # Imports the solver RADAU from Assimulo G_EARTH = np.array([0.0, 0.0, -9.81]) # gravitational acceleration # Example one: # Falling mass. # State vector y = mass.pos, mass.vel # Derivative yd = mass.vel, mass.acc # Residual res = (y.vel - yd.vel), (yd.acc - G_EARTH) def res1(t, y, yd): res_0 = y[3:6] - yd[0:3] res_1 = yd[3:6] - G_EARTH return np.append(res_0, res_1) def run_example(): # Set the initial conditons t0 = 0.0 # Initial time vel_0 = [0.0, 0.0, 50.0] # Initial velocity pos_0 = [0.0, 0.0, 0.0] # Initial position acc_0 = [0.0, 0.0, -9.81] # Initial acceleration y0 = pos_0 + vel_0 # Initial pos, vel yd0 = vel_0 + acc_0 # Initial vel, acc model = Implicit_Problem(res1, y0, yd0, t0) # Create an Assimulo problem model.name = 'Falling mass' # Specifies the name of problem (optional) sim = Radau5DAE(model) # Create the Radau solver tfinal = 10.0 # Specify the final time ncp = 500 # Number of communcation points (number of return points) # Use the .simulate method to simulate and provide the final time and ncp (optional) time, y, yd = sim.simulate(tfinal, ncp) # plot the result pos_z = y[:,2] vel_z = y[:,5]
The Python package Assimulo has a very good documentation for solving implicit problems:
I did not found anything similar in the documentation of DifferentialEquations.
Which problem shall I define to solve this example?
Shall I create a DAEProblem? But then I cannot use Radau, correct?
Please keep in mind that my real problem is large, I coded it in Numba/Python (now converting it to Julia) and would prefer to stick to the Radau solver, because it performs well for my stiff problem of about 50 equations.