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:
https://jmodelica.org/assimulo/tutorial_imp.html
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.
