I want to introduce Julia to one course. We have been using MATLAB in this course for many years, which are very robust. For these special problems, I tried Python and it seems not as robust as MATLAB (I am also a beginner of Python), but I managed to solve them all. I heard from friends about Julia, and I tried, but apparently I can hardly solve the same problems (the warning/error is “Instability detected. Aborting”). Here below is the source code, and it seems only the solver
Rosenbrock23 gives the expected results. This situation makes it difficult for me to explain to others (in MATLAB almost all solvers work, just the efficiency matters).
using DifferentialEquations, Plots function dyth_dtau(dyth, yth, par, t) # nt = int(length(yth)/2) nt = par # this can be calculated from the line given above dz = par # this can be calculated from n (nt) gamma = par beta = par phi = par Le = par theta = yth[1:nt] y = yth[nt+1:2*nt] rhs = [phi^2 * y[i] * exp(gamma*(1.0 - 1.0/theta[i])) for i in 1:nt] # dtheta/dtau dyth = 2 * (theta - theta)/dz^2 + beta * rhs for i in 2:nt-1 dyth[i] = (theta[i+1] - 2 * theta[i] + theta[i-1])/dz^2 + beta * rhs[i] end dyth[nt] = (1 - 2 * theta[nt] + theta[nt-1])/dz^2 + beta * rhs[nt] dyth[1:nt] = dyth[1:nt] / Le # dy/dtau dyth[nt+1] = 2 * (y - y)/dz^2 - rhs for i in 2:nt-1 dyth[nt+i] = (y[i+1] - 2 * y[i] + y[i-1])/dz^2 - rhs[i] end dyth[2*nt] = (1 - 2 * y[nt] + y[nt-1])/dz^2 - rhs[nt] end n = 99 dz = 1/(n+1) gamma = 30.0 beta = 0.15 phi = 1.1 Le_arr = [1.0, 0.5, 0.25, 0.15, 0.1] t_span = (0.0, 1.0) yth0 = ones(2*n) for i in 1:length(Le_arr) Le = Le_arr[i] println(Le) par = (n, dz, gamma, beta, phi, Le) prob = ODEProblem(dyth_dtau, yth0, t_span, par) sol = solve(prob) if (i == 1) global plt=plot(sol.t, sol[1, :]) else global plt=plot!(sol.t, sol[1, :]) end end display(plt)
Here is the status of the package: " [0c46a032] DifferentialEquations v6.16.0 ".
It is very strange that not all ODE solvers are available. For example: CVODE_BDF() and RadauIIA() are missing.
Could you please give me some advices? Thank you.