How to convert the following mathematica code to Matlab?

3 visualizzazioni (ultimi 30 giorni)
I am inexperienced in Matlab and need to convert the following mathematica codes in it. Any help will be appreciated.
ClearAll["Global`*"];
tstar = -5;
ode1 = y'[t] == -Sin[x[t]]/y[t];
ode2 = x'[ t] == -Cos[x[t]] (6 Sin[x[t]] Cos[x[t]] + y[t] (b - c (1 + 3*y[t]^2)))/(2*y[t]^3*(b + c (y[t]^2 - 1)));
ode3 = v'[t] == -(b + c*(y[t]^2 - 1))/(4*y[t]*Cos[x[t]]) + Sin[x[t]]/(2*y[t]^2);
bc = {x[tstar] == 0, y[tstar] == Br, v[tstar] == Log[Dr]};
sol = ParametricNDSolve[{ode1, ode2, ode3, bc}, {x, y, v}, {t, tstar, 0}, {b, c, Br, Dr}]
data = Table[{Br, Dr} /. FindRoot[{(y[b, c, Br, Dr][0] - 1) /. sol, v[b, c, Br, Dr][0] /. sol}, {{b, 1}, {c, 1}}, WorkingPrecision -> 5], {Dr, 20, 50, 1}, {Br, 1.5, 2.7, .1}]
  2 Commenti
Fabio Freschi
Fabio Freschi il 18 Ott 2022
The best thing to do is to understand the math and programming natively in the new language. I am not a Mathematica user, so I cannot help unless you find out a way to describe your problem with "math" language, and not with mathematica syntax
Fun123
Fun123 il 18 Ott 2022
Modificato: Fun123 il 18 Ott 2022
Than you so much. I have three odes(ode1,ode2,ode3) with four parameters( Br, Dr, b and c). Mathematica code parametricndsolve provides parametric solution for the initail conditions(here named as bc) at tstar=-5. Then in the last line (data =Table....) we run two loops for parameters Dr=20 to 50 with stepsize 1 and then Br=1.5 to 2.7 with step size 0.1 to determine other two parameters b and c by using mathematica code Findroot for the boundary condition v(0)=0 and y(0)=1.

Accedi per commentare.

Risposte (1)

Sam Chak
Sam Chak il 18 Ott 2022
Identify the type of problem at first. If the problem is non-stiff, use ode45. Miscellaneous examples are provided in the links.
help ode45
ODE45 Solve non-stiff differential equations, medium order method. [TOUT,YOUT] = ODE45(ODEFUN,TSPAN,Y0) integrates the system of differential equations y' = f(t,y) from time TSPAN(1) to TSPAN(end) with initial conditions Y0. Each row in the solution array YOUT corresponds to a time in the column vector TOUT. * ODEFUN is a function handle. For a scalar T and a vector Y, ODEFUN(T,Y) must return a column vector corresponding to f(t,y). * TSPAN is a two-element vector [T0 TFINAL] or a vector with several time points [T0 T1 ... TFINAL]. If you specify more than two time points, ODE45 returns interpolated solutions at the requested times. * YO is a column vector of initial conditions, one for each equation. [TOUT,YOUT] = ODE45(ODEFUN,TSPAN,Y0,OPTIONS) specifies integration option values in the fields of a structure, OPTIONS. Create the options structure with odeset. [TOUT,YOUT,TE,YE,IE] = ODE45(ODEFUN,TSPAN,Y0,OPTIONS) produces additional outputs for events. An event occurs when a specified function of T and Y is equal to zero. See ODE Event Location for details. SOL = ODE45(...) returns a solution structure instead of numeric vectors. Use SOL as an input to DEVAL to evaluate the solution at specific points. Use it as an input to ODEXTEND to extend the integration interval. ODE45 can solve problems M(t,y)*y' = f(t,y) with mass matrix M that is nonsingular. Use ODESET to set the 'Mass' property to a function handle or the value of the mass matrix. ODE15S and ODE23T can solve problems with singular mass matrices. ODE23, ODE45, ODE78, and ODE89 are all single-step solvers that use explicit Runge-Kutta formulas of different orders to estimate the error in each step. * ODE45 is for general use. * ODE23 is useful for moderately stiff problems. * ODE78 and ODE89 may be more efficient than ODE45 on non-stiff problems that are smooth except possibly for a few isolated discontinuities. * ODE89 may be more efficient than ODE78 on very smooth problems, when integrating over long time intervals, or when tolerances are tight. Example [t,y]=ode45(@vdp1,[0 20],[2 0]); plot(t,y(:,1)); solves the system y' = vdp1(t,y), using the default relative error tolerance 1e-3 and the default absolute tolerance of 1e-6 for each component, and plots the first component of the solution. Class support for inputs TSPAN, Y0, and the result of ODEFUN(T,Y): float: double, single See also ODE23, ODE78, ODE89, ODE113, ODE15S, ODE23S, ODE23T, ODE23TB, ODE15I, ODESET, ODEPLOT, ODEPHAS2, ODEPHAS3, ODEPRINT, DEVAL, ODEEXAMPLES, FUNCTION_HANDLE. Documentation for ode45 doc ode45 Other uses of ode45 dlarray/ode45
Try ode15s when ode45 fails or is inefficient and you suspect that the problem is stiff.
help ode15s
ODE15S Solve stiff differential equations and DAEs, variable order method. [TOUT,YOUT] = ODE15S(ODEFUN,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates the system of differential equations y' = f(t,y) from time T0 to TFINAL with initial conditions Y0. ODEFUN is a function handle. For a scalar T and a vector Y, ODEFUN(T,Y) must return a column vector corresponding to f(t,y). Each row in the solution array YOUT corresponds to a time returned in the column vector TOUT. To obtain solutions at specific times T0,T1,...,TFINAL (all increasing or all decreasing), use TSPAN = [T0 T1 ... TFINAL]. [TOUT,YOUT] = ODE15S(ODEFUN,TSPAN,Y0,OPTIONS) solves as above with default integration properties replaced by values in OPTIONS, an argument created with the ODESET function. See ODESET for details. Commonly used options are scalar relative error tolerance 'RelTol' (1e-3 by default) and vector of absolute error tolerances 'AbsTol' (all components 1e-6 by default). If certain components of the solution must be non-negative, use ODESET to set the 'NonNegative' property to the indices of these components. The 'NonNegative' property is ignored for problems where there is a mass matrix. The Jacobian matrix df/dy is critical to reliability and efficiency. Use ODESET to set 'Jacobian' to a function handle FJAC if FJAC(T,Y) returns the Jacobian df/dy or to the matrix df/dy if the Jacobian is constant. If the 'Jacobian' option is not set (the default), df/dy is approximated by finite differences. Set 'Vectorized' 'on' if the ODE function is coded so that ODEFUN(T,[Y1 Y2 ...]) returns [ODEFUN(T,Y1) ODEFUN(T,Y2) ...]. If df/dy is a sparse matrix, set 'JPattern' to the sparsity pattern of df/dy, i.e., a sparse matrix S with S(i,j) = 1 if component i of f(t,y) depends on component j of y, and 0 otherwise. ODE15S can solve problems M(t,y)*y' = f(t,y) with mass matrix M(t,y). Use ODESET to set the 'Mass' property to a function handle MASS if MASS(T,Y) returns the value of the mass matrix. If the mass matrix is constant, the matrix can be used as the value of the 'Mass' option. Problems with state-dependent mass matrices are more difficult. If the mass matrix does not depend on the state variable Y and the function MASS is to be called with one input argument T, set 'MStateDependence' to 'none'. If the mass matrix depends weakly on Y, set 'MStateDependence' to 'weak' (the default) and otherwise, to 'strong'. In either case the function MASS is to be called with the two arguments (T,Y). If there are many differential equations, it is important to exploit sparsity: Return a sparse M(t,y). Either supply the sparsity pattern of df/dy using the 'JPattern' property or a sparse df/dy using the Jacobian property. For strongly state-dependent M(t,y), set 'MvPattern' to a sparse matrix S with S(i,j) = 1 if for any k, the (i,k) component of M(t,y) depends on component j of y, and 0 otherwise. If the mass matrix is non-singular, the solution of the problem is straightforward. See examples FEM1ODE, FEM2ODE, BATONODE, or BURGERSODE. If M(t0,y0) is singular, the problem is a differential- algebraic equation (DAE). ODE15S solves DAEs of index 1. DAEs have solutions only when y0 is consistent, i.e., there is a yp0 such that M(t0,y0)*yp0 = f(t0,y0). Use ODESET to set 'MassSingular' to 'yes', 'no', or 'maybe'. The default of 'maybe' causes ODE15S to test whether M(t0,y0) is singular. You can provide yp0 as the value of the 'InitialSlope' property. The default is the zero vector. If y0 and yp0 are not consistent, ODE15S treats them as guesses, tries to compute consistent values close to the guesses, and then goes on to solve the problem. See examples HB1DAE or AMP1DAE. [TOUT,YOUT,TE,YE,IE] = ODE15S(ODEFUN,TSPAN,Y0,OPTIONS) with the 'Events' property in OPTIONS set to a function handle EVENTS, solves as above while also finding where functions of (T,Y), called event functions, are zero. For each function you specify whether the integration is to terminate at a zero and whether the direction of the zero crossing matters. These are the three column vectors returned by EVENTS: [VALUE,ISTERMINAL,DIRECTION] = EVENTS(T,Y). For the I-th event function: VALUE(I) is the value of the function, ISTERMINAL(I)=1 if the integration is to terminate at a zero of this event function and 0 otherwise. DIRECTION(I)=0 if all zeros are to be computed (the default), +1 if only zeros where the event function is increasing, and -1 if only zeros where the event function is decreasing. Output TE is a column vector of times at which events occur. Rows of YE are the corresponding solutions, and indices in vector IE specify which event occurred. SOL = ODE15S(ODEFUN,[T0 TFINAL],Y0...) returns a structure that can be used with DEVAL to evaluate the solution or its first derivative at any point between T0 and TFINAL. The steps chosen by ODE15S are returned in a row vector SOL.x. For each I, the column SOL.y(:,I) contains the solution at SOL.x(I). If events were detected, SOL.xe is a row vector of points at which events occurred. Columns of SOL.ye are the corresponding solutions, and indices in vector SOL.ie specify which event occurred. Example [t,y]=ode15s(@vdp1000,[0 3000],[2 0]); plot(t,y(:,1)); solves the system y' = vdp1000(t,y), using the default relative error tolerance 1e-3 and the default absolute tolerance of 1e-6 for each component, and plots the first component of the solution. See also ODE23S, ODE23T, ODE23TB, ODE45, ODE23, ODE113, ODE15I, ODESET, ODEPLOT, ODEPHAS2, ODEPHAS3, ODEPRINT, DEVAL, ODEEXAMPLES, VDPODE, BRUSSODE, HB1DAE, FUNCTION_HANDLE. Documentation for ode15s doc ode15s
If it is a boundary value problem, then use bvp4c.
help bvp4c
BVP4C Solve boundary value problems for ODEs by collocation. SOL = BVP4C(ODEFUN,BCFUN,SOLINIT) integrates a system of ordinary differential equations of the form y' = f(x,y) on the interval [a,b], subject to general two-point boundary conditions of the form bc(y(a),y(b)) = 0. ODEFUN and BCFUN are function handles. For a scalar X and a column vector Y, ODEFUN(X,Y) must return a column vector representing f(x,y). For column vectors YA and YB, BCFUN(YA,YB) must return a column vector representing bc(y(a),y(b)). SOLINIT is a structure with fields x -- ordered nodes of the initial mesh with SOLINIT.x(1) = a, SOLINIT.x(end) = b y -- initial guess for the solution with SOLINIT.y(:,i) a guess for y(x(i)), the solution at the node SOLINIT.x(i) BVP4C produces a solution that is continuous on [a,b] and has a continuous first derivative there. The solution is evaluated at points XINT using the output SOL of BVP4C and the function DEVAL: YINT = DEVAL(SOL,XINT). The output SOL is a structure with SOL.solver -- 'bvp4c' SOL.x -- mesh selected by BVP4C SOL.y -- approximation to y(x) at the mesh points of SOL.x SOL.yp -- approximation to y'(x) at the mesh points of SOL.x SOL.stats -- computational cost statistics (also displayed when the 'Stats' option is set with BVPSET). SOL = BVP4C(ODEFUN,BCFUN,SOLINIT,OPTIONS) solves as above with default parameters replaced by values in OPTIONS, a structure created with the BVPSET function. To reduce the run time greatly, use OPTIONS to supply a function for evaluating the Jacobian and/or vectorize ODEFUN. See BVPSET for details and SHOCKBVP for an example that does both. Some boundary value problems involve a vector of unknown parameters p that must be computed along with y(x): y' = f(x,y,p) 0 = bc(y(a),y(b),p) For such problems the field SOLINIT.parameters is used to provide a guess for the unknown parameters. On output the parameters found are returned in the field SOL.parameters. The solution SOL of a problem with one set of parameter values can be used as SOLINIT for another set. Difficult BVPs may be solved by continuation: start with parameter values for which you can get a solution, and use it as a guess for the solution of a problem with parameters closer to the ones you want. Repeat until you solve the BVP for the parameters you want. The function BVPINIT forms the guess structure in the most common situations: SOLINIT = BVPINIT(X,YINIT) forms the guess for an initial mesh X as described for SOLINIT.x, and YINIT either a constant vector guess for the solution or a function handle. If YINIT is a function handle then for a scalar X, YINIT(X) must return a column vector, a guess for the solution at point x in [a,b]. If the problem involves unknown parameters SOLINIT = BVPINIT(X,YINIT,PARAMS) forms the guess with the vector PARAMS of guesses for the unknown parameters. BVP4C solves a class of singular BVPs, including problems with unknown parameters p, of the form y' = S*y/x + f(x,y,p) 0 = bc(y(0),y(b),p) The interval is required to be [0, b] with b > 0. Often such problems arise when computing a smooth solution of ODEs that result from PDEs because of cylindrical or spherical symmetry. For singular problems the (constant) matrix S is specified as the value of the 'SingularTerm' option of BVPSET, and ODEFUN evaluates only f(x,y,p). The boundary conditions must be consistent with the necessary condition S*y(0) = 0 and the initial guess should satisfy this condition. BVP4C can solve multipoint boundary value problems. For such problems there are boundary conditions at points in [a,b]. Generally these points represent interfaces and provide a natural division of [a,b] into regions. BVP4C enumerates the regions from left to right (from a to b), with indices starting from 1. In region k, BVP4C evaluates the derivative as YP = ODEFUN(X,Y,K). In the boundary conditions function, BCFUN(YLEFT,YRIGHT), YLEFT(:,K) is the solution at the 'left' boundary of region k and similarly for YRIGHT(:,K). When an initial guess is created with BVPINIT(XINIT,YINIT), XINIT must have double entries for each interface point. If YINIT is a function handle, BVPINIT calls Y = YINIT(X,K) to get an initial guess for the solution at X in region k. In the solution structure SOL returned by BVP4C, SOL.x has double entries for each interface point. The corresponding columns of SOL.y contain the 'left' and 'right' solution at the interface, respectively. See THREEBVP for an example of solving a three-point BVP. Example solinit = bvpinit([0 1 2 3 4],[1 0]); sol = bvp4c(@twoode,@twobc,solinit); solve a BVP on the interval [0,4] with differential equations and boundary conditions computed by functions twoode and twobc, respectively. This example uses [0 1 2 3 4] as an initial mesh, and [1 0] as an initial approximation of the solution components at the mesh points. xint = linspace(0,4); yint = deval(sol,xint); evaluate the solution at 100 equally spaced points in [0 4]. The first component of the solution is then plotted with plot(xint,yint(1,:)); For more examples see TWOBVP, FSBVP, SHOCKBVP, MAT4BVP, EMDENBVP, THREEBVP. See also BVP5C, BVPSET, BVPGET, BVPINIT, BVPXTEND, DEVAL, FUNCTION_HANDLE. Documentation for bvp4c doc bvp4c

Community Treasure Hunt

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

Start Hunting!

Translated by