Translating from python to matlab code

5 views (last 30 days)
Alex Monserrat
Alex Monserrat on 29 Apr 2022
Commented: Image Analyst on 29 Apr 2022
I want to translate this python code which solves an ODE with radau mode to matlab. I am new to matlab. There is no documentation about the radau mode, can somebody give me a hand ?
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
def deriv(t, y):
"""ODEs for Robertson's chemical reaction system."""
x, y, z = y
xdot = -0.04 * x + 1.e4 * y * z
ydot = 0.04 * x - 1.e4 * y * z - 3.e7 * y**2
zdot = 3.e7 * y**2
return xdot, ydot, zdot
# Initial and final times.
t0, tf = 0, 500
# Initial conditions: [X] = 1; [Y] = [Z] = 0. # Given as a TUPLE
y0 = 1, 0, 0
# Solve, using a method resilient to stiff ODEs.
soln = solve_ivp(deriv, (t0, tf), y0, method='Radau')
print(soln.nfev, 'evaluations required.')
# print(soln) all results
# Plot the concentrations as a function of time. Scale [Y] by 10**YFAC
# so its variation is visible on the same axis used for [X] and [Z].
YFAC = 4
plt.plot(soln.t, soln.y[0], label='[X]')
plt.plot(soln.t, 10**YFAC*soln.y[1], label=r'$10^{}\times$[Y]'.format(YFAC))
plt.plot(soln.t, soln.y[2], label='[Z]')
plt.xlabel('time /s')
plt.ylabel('concentration /arb. units')
plt.legend()
plt.show()

Answers (1)

Image Analyst
Image Analyst on 29 Apr 2022
Replace # with %.
Get rid of the plt., for example plt.xlabel goes to simply xlabel.
No need for the first 3 lines of code (import and from).
printf() goes to fprintf(), like fprintf('%f\n', yourNumber). Look up the documentation.
plt.plot(soln.t, soln.y[0], label='[X]')
would be like
plot(soln.t, soln.y, 'b-');
No need for the show() function.
Try with that then tackle the remaining problems one at a time.
  2 Comments
Alex Monserrat
Alex Monserrat on 29 Apr 2022
Ok, thank you.
How do I mount the ode system and solver ?
function radau_try
tspan = [0 500];
x0 = 1;
y0 = 0;
z0 = 0;
dx(1) = -0.04 * x + 1.e4 * y * z;
dy(1) = 0.04 * x - 1.e4 * y * z - 3.e7 * y^2;
dz(1) = 3.e7 * y^2;
% [t,y] = ode113(@(t,y) 2*t, tspan, y0);
radausolver(OdeFcn,tspan,y,Op,varargin)
end
Image Analyst
Image Analyst on 29 Apr 2022
Sorry, I don't know -- I don't use those functions.

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!