Translating from python to matlab code

5 visualizzazioni (ultimi 30 giorni)
Alex Monserrat
Alex Monserrat il 29 Apr 2022
Commentato: Image Analyst il 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()

Risposte (1)

Image Analyst
Image Analyst il 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 Commenti
Alex Monserrat
Alex Monserrat il 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 il 29 Apr 2022
Sorry, I don't know -- I don't use those functions.

Accedi per commentare.

Categorie

Scopri di più su Programming in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by