Translating from python to matlab code

4 vues (au cours des 30 derniers jours)
Alex Monserrat
Alex Monserrat le 29 Avr 2022
Commenté : Image Analyst le 29 Avr 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()

Réponses (1)

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

Connectez-vous pour commenter.

Catégories

En savoir plus sur Programming dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by