Tag: differential

Getting Mathematica to solve a system of two second order nonlinear ordinary differential equations

I tried solving a system of two second order nonlinear ordinary differential equations using the DSolve command. First, I tried like this:

eqns = {A''[x] == 2/B[x]*A'[x]*B'[x],     B''[x] + 1/B[x]*(A'[x])^2 - 1/B[x]*(B'[x])^2 == 0}; sol = DSolve[eqns, {A, B}, x]

However, as Mathematica didn’t (couldn’t?) solve this, I transformed it into a system of four first order equations:

eqns = {c'[x] == 2/B[x]*c[x]*d[x],     d'[x] + 1/B[x]*(c[x])^2 - 1/B[x]*(d[x])^2 == 0, c[x] == A'[x],     d[x] == B'[x], c[0] == 1, d[0] == 1, A[0] == 1, B[0] == 1}; sol = DSolve[eqns, {A, B, c, d}, x]

This still doesn’t work. Weirdly enough, I don’t even get an error message.

I only included the boundary conditions thinking that they may be helpful, but they aren’t part of my original problem.

Your help would be greatly appreciated:)

Fit data to a model of two-variable differential equations

I am trying to fit some data to the solution of a differential equation. I am using "ParametricNDSolve" to create a model which should fit the data using "FindFit".

Now, the code I propose is very similar to the examples I have seen to solve similar problems. However, the equation of this case is a little more complicated that the examples of Mathematica manuals (because is a function of two variables, and the initials conditions involves two numeric integrals).

The code keeps running for hours, and it does not appear to solve the problems at all:

data = {{21/100, 0.260882276365091`}, {6/25, 0.330910580727009`}, {27/ 100, 0.271601829283246`}, {3/10, 0.294039749043066`}, {33/100,  0.373363616981994`}, {9/25, 0.467495450916`}, {39/100,  0.50972503848512`}, {21/50, 0.639300915026114`}, {9/20,  0.679672806174314`}, {12/25, 0.693446859556429`}, {51/100,  0.697043731283207`}, {27/50, 0.736147748563448`}, {57/100,  0.706580484286456`}, {3/5, 0.719128578284264`}, {63/100,  0.762276173639189`}, {33/50, 0.788753648395851`}, {69/100,  0.802107836002146`}, {18/25, 0.803992011056852`}, {3/4,  0.796349018322968`}, {39/50, 0.795845629380594`}, {81/100,  0.808805024532776`}, {21/25, 0.812210415525446`}, {87/100,  0.811677706096654`}, {9/10, 0.796614634731033`}, {93/100,  0.802573366818874`}, {24/25, 0.801851431134633`}, {99/100,  0.799762250401537`}, {51/50, 0.803234757269179`}, {21/20,  0.812346679044962`}, {27/25, 0.804974524275207`}, {111/100,  0.813054479870107`}, {57/50, 0.802494974488922`}, {117/100,  0.811124462792484`}, {6/5, 0.810259975308694`}, {123/100,  0.810910967136978`}, {63/50, 0.815891670404585`}, {33/25,  0.803843274767398`}, {69/50, 0.814778081357444`}, {36/25,  0.809686831953762`}, {3/2, 0.829532613157218`}, {9/5,  0.815351190536115`}, {21/10, 0.828143435644695`}, {12/5,  0.791269047273677`}, {27/10, 0.821321230878138`}, {3,  0.81153212440728`}, {33/10, 0.821798500231733`}, {18/5,  0.819594302125574`}, {39/10, 0.828629593518031`}};  model = ParametricNDSolveValue[{ D[v[x, t], t] == d1*D[v[x, t], x, x], v[-5, t] == 0, v[10, t] == flux, v[x, 0] == flux*NIntegrate[Cos[y]^n1, {y, ArcTan[-b1*(x + a1)], Pi/2}]/NIntegrate[Cos[y]^n1, {y, -Pi/2, Pi/2}]}, v, {x, -5, 10}, {t, 0, 2}, {d1, flux, n1, b1, a1}]  fit = FindFit[data,  (c*model[d1, flux, n1, b1, a1][x, 2])/((c - 1)*model[d1, flux, n1, b1, a1][x, 2] + 1), {{d1, 0.001}, {flux, 0.57}, {c, 3.2}, {n1, 2}, {b1,1}, {a1, -0.4}}, x]

I guess that this code is computationally inefficient, but I need help to understand why.

Thank you all very much

Ps: english is obviously not my first lenguage, so I apologize for the grammar mistakes.

How to minimize a problem with differential equations in it?

My problem is to minimize a target function with differential equations in it. The target function is $ $ \min \mathrm{J}\left( p \right) =\sum_{i=1}^5{|x_i\left( p \right) -y_i|}\,\,$ $ where $ x_i(p)$ are differential equations like this:

stateFunctionV[v1_, v2_, v3_, v4_, v5_, v6_, v7_, v8_] := { \[Mu] = (67/100)*(x2[t]/(x2[t] + 28/100))*(1 - x3[t]/940)*(1 - x4[t]/1050)^3*(1 - x5[t]/361); x1'[t] == \[Mu]*x1[t], x2'[t] == -(v1 + \[Mu]/v5)* x1[t], x3'[t] == (v2 + \[Mu]*v6)*x1[t], x4'[t] == (v3 + \[Mu]*v7)* x1[t], x5'[t] == (v4 + \[Mu]*v8)* x1[t], x1[0] == 0.2025, x2[0] == 441.337, x3[0] == 0, x4[0] == 0,  x5[0] == 0}

and $ y_i$ is the experiment data.

The code:

PIP1[value_] := { Sum[Sum[Abs[ NDSolveValue[{ [Mu] = (67/100)*(x2[t]/(x2[t] + 28/100))*(1 - x3[t]/940)*(1 - x4[t]/1050)^3*(1 - x5[t]/361); x1'[t] == \[Mu]*x1[t], x2'[t] == -(value[[1]] + \[Mu]/value[[5]])* x1[t], x3'[t] == (value[[2]] + \[Mu]*value[[6]])*x1[t], x4'[t] == (value[[3]] + \[Mu]*value[[7]])* x1[t], x5'[t] == (value[[4]] + \[Mu]*value[[8]])* x1[t], x1[0] == 0.2025, x2[0] == 441.337, x3[0] == 0, x4[0] == 0,  x5[0] == 0}, {x1, x2, x3, x4, x5}, {t, 0, 6.5}][[i]] [data2[[i, j, 1]]] - data2[[i, j, 2]] ], {j, 1, 7}], {i, 1, 5}]} NMinimize[PIP1[{v1, v2, v3, v4, v5, v6, v7, v8}], {v1, v2, v3, v4, v5, v6, v7, v8}]

doesn’t work. I also tried to take the differential functions out as a new function and it also failed.

So I wonder if there’s a way to minimize this problem with NMinimize or not? I would be grateful for your advice on how to couple NMinimize and NDSolveValue.

The experiment data are below:

data2 = {{{0, 0.2025}, {2, 0.76}, {3, 1.28}, {4, 2.2}, {5, 2.75}, {6,       2.95}, {6.5, 3.055}}, {{0, 441.337}, {2, 400.435}, {3,       300.326}, {4, 219.348}, {5, 120.109}, {6, 50}, {6.5, 25}}, {{0,       0}, {2, 42.2895}, {3, 78.6579}, {4, 117.026}, {5, 152.158}, {6,       229.553}, {6.5, 224.1205}}, {{0, 0}, {2, 1.66667}, {3,       13.55}, {4, 16.6667}, {5, 22.7333}, {6, 34.4333}, {6.5,       71.2333}}, {{0, 0}, {2, 6.73913}, {3, 7.86957}, {4, 9.3913}, {5,       8.34783}, {6, 10.9565}, {6.5, 11.8261}}};

Full line trajectories plot for the solution of Second Order nonlinear coupled differential equations

I wanted to plot a phase plane containing the trajectories of the solutions found by using ‘NDSolve’ using the initial conditions for x[0], y[0], x'[0] and y'[0]. The equations are: x”[t] – 2 y'[t] == -x[t] + y[t]^2; y”[t] + 2 x'[t] == x[t] + y[t] + x[t]*y[t]

The equilibrium point for the system is (0,0). I have plotted the stream plot for the system but unable to plot a phase portrait that would give me the full line trajectories of the system for different initial conditions. I am also looking for any periodic solution if present in it. The stream plot I got is given below and I would take initial conditions from it.

enter image description here

I get this by using the Parametric Plot of the NDSolve solution:

enter image description here

Kindly help in this capacity. Thanks in advance.

Nonlinear Differential Equation with Paramaters in Sqrt Function

For a university project, I am trying to see if my system will have choked flow and also plot the resulting pressure spike. I set up the system below to try to model the transient response.

I am able to solve the equations TEq and mEq simultaneously if I use the mDotOutChoked equation in place of the mDotOut equation. But I want to know if it ever even reaches the choked condition, and what the steady-state pressure and temperature will be so I want to start with just using the mDotOut equation.

When I run this as it is below, I get "This computation has exceeded the time limit for your plan" followed by "{sol1} is neither a list or replacement rules nor a valid dispatch table, and so…" for all the graphs.

If I replace the pressure term with just a constant, I can get some limited success depending on the value I use. This makes me think there might be an issue with having the T[t] and P[t] terms in the square root in the mDotOut equation.

Is there some issue in {sol1} or is that issue just a result of me needing to upgrade for more compute power?

ClearAll["Global`*"]  Q =  100 ;(*Heat into the shroud in Watts. Based on roughly 1350 W/m^2 from the solar simulator on one face of the shroud*) QAmb = 0  ;(*Heat loss to ambinet. Zero for now*) A =  2*1.935*10^-5 ;(*Area of orifice m^2  Based on 1/4 inch pipe with 0.028 inch wall thickness (two outlets)*) h =  199000 ;(*Heat of vaporization of LN2 J/kg*) Cp =  1039 ;(*Specific heat of nitrogen J/(kg K) *) R =  296.8 ;(*Gas constant for nitrogen J/(kg K)*) γ = 1.40  ;(*Specific heat ratio*) V = 0.001 ; (*Enclosed volume m^3*) Pe = 101000 ; (*External pressure in pa*) ρo = 4 ;(*Approx density of nitrogen at 80K in kg/m^3. This was the lowest temp data I could find*) tf = 300 ;(*Final time in seconds*)  P = m[t]*R*T[t]/V ;(*Pressure term*) mDotEvap = Q/h ; (*rate of evap*) mDotOut = (P*A/Sqrt[T[t]])*Sqrt[(2γ/(R(γ-1)))*((Pe/P)^(2/γ)-(Pe/P)^((γ+1)/γ))];  (*mass flow out of the orifice*) mDotOutChoked =  (P*A/Sqrt[T[t]])*Sqrt[γ/R]*(2/(γ+1))^((γ+1)/(2(γ-1))); (*mass flow out of the orifice if choked*)  TEq = T'[t] == 1/(m[t]*Cp)(mDotEvap*h-mDotOut*Cp*T[t]-QAmb) ; (*Diff Eq for Temperature in the cavity*) mEq = m'[t] == mDotEvap - mDotOut ; (*Conservation of mass*)  icT = T[0] == 77 ;(*initial temp in the cavity in K*) icm = m[0] == ρo*V ;(*initial mass of the vaporized gas. Assuming it just starts at 77k at 1atm and then adding heat*)  sol1 = NDSolve[{TEq,mEq, icT, icm}, {T[t], m[t]},{t, 0, tf} ] ; P2[t_] = m[t]*R*T[t]/V /.sol1 ; (*Plugging back to get shroud pressure as functon of time*)  Plot[{T[t]/.sol1},{t,0,tf},PlotRange -> Automatic, ImageSize->"Large",PlotLabels->Automatic, AxesLabel -> {Time (s),Temperature (K) }] Plot[{m[t]/.sol1},{t,0,tf},PlotRange -> Automatic, ImageSize->"Large",PlotLabels->Automatic,  AxesLabel -> {Time (s),Mass (Kg) }] Plot[P2[t], {t, 0, tf}, PlotRange -> Automatic, ImageSize -> "Large", PlotLabels -> Automatic,  AxesLabel -> {Time (s),Pressure (Pa) }]

NMinimize differential evolution: how do search points and initial points really work?

I am puzzled by the settings of a specific method for NMinimize called DifferentialEvolution, given for example here. The settings include both "SearchPoints", or the size of the population of evolving points, and "InitialPoints", or the initial population(?). What I don’t quite understand is that "SearchPoints" can be set to a value different from the number of specified "InitialPoints" and NMinimize will often happily proceed, and sometimes it won’t.

What might be happening under the hood? If there are less "InitialPoints" than there are "SearchPoints", will the population be filled up with random points to meet the specified number of "SearchPoints" before evolution begins or does it proceed with a population that is in fact smaller than "SearchPoints"? What about the opposite case?

Here is an example where we give much less "InitialPoints" than there are "SearchPoints".

Clear[f, c, v, x1, x2, y1, y2, y3]; f = 2 x1 + 3 x2 + 3 y1/2 + 2 y2 - y3/2; c = {x1^2 + y1 == 5/4, x2^(3/2) + 3 y2/2 == 3, x1 + y1 <= 8/5,     4 x2/3 + y2 <= 3, y3 <= y1 + y2, 0 <= x1 <= 10, 0 <= x2 <= 10,     0 <= y1 <= 1, 0 <= y2 <= 1,     0 <= y3 <= 1, {y1, y2, y3} \[Element] Integers}; v = {x1, x2, y1, y2, y3};  NMinimize[{f, c}, v, Method -> "DifferentialEvolution"] (*{7.66718, {x1 -> 1.11803, x2 -> 1.31037, y1 -> 0, y2 -> 1, y3 -> 1}}*)  points = 5; searchpoints = 50; listpoints = {RandomReal[{0, 10}, points],  RandomReal[{0, 10}, points], RandomReal[{0, 1}, points],  RandomReal[{0, 1}, points],  RandomReal[{0, 1}, points]}\[Transpose];  NMinimize[{f, c}, v,   Method -> {"DifferentialEvolution", "SearchPoints" -> searchpoints,  "InitialPoints" -> listpoints}] (*{7.66718, {x1 -> 1.11803, x2 -> 1.31037, y1 -> 0, y2 -> 1, y3 -> 1}}*)

If there are more "InitialPoints" than there are "SearchPoints", NMinimize sometimes works and sometimes doesn’t, depending on the value of the random seed for instance.

Can built-in functions deal with stochastic delay differential equations (SDDE)?

I know that functions like NDSolve can deal with delay differential equations and in the meanwhile, functions like ItoProcess and RandomFunction handle stochastic differential equations. So I wonder whether any built-in functions can handle it when the above two cases are combined together. For example, I naively tried the below codes by just slightly modifying the first example of ItoProcess (x[t] -> x[t - 1] in the square root)

proc = ItoProcess[\[DifferentialD]x[t] == -x[t] \[DifferentialD]t + Sqrt[1 + x[t - 1]^2]         \[DifferentialD]w[t], x[t], {x, 1}, t, w \[Distributed] WienerProcess[]] RandomFunction[proc, {0., 5., 0.01}]

The first row of codes runs seemly well, but the second one just returns a RandomFunction::unsproc error, specifically RandomFunction::unsproc: The specification `<Ito process>` is not a random process recognized by the system..

Or do I have to implement a version myself with Euler method alike?

How do i use Streamplot to plot a non homogenous differential equation

I have the following equations:

x’ = x-y y’ = x+y-2xy

I used the following code to do the streamplot: I apologize for only having an image, as I am using Mathematica through my school servers I cannot copy and paste the data.

enter image description here

But I do not get any results, so then I linearize the equations and plot them separately, but I cannot figure out how to combine the plots to show in one plot.:

enter image description here

Is there any way to use Streamplot with nonhomogenous DE’s or a better way to combine the equations?

Why both DSolve and NDSolve are unable to solve a second-order differential equation?

I am trying to solve a recurrence relation using generating functions method. After some long calculations, I have arrived to this second-order differential equation: $ 0.5 x^5 y”(x)+(2x^4+x^3)y'(x)+\left(x^3+x^2+x-1\right)y(x)+1=0$

and these conditions: $ y(0)=1, y'(0)=1$ . $ y(x)$ is the function that needs to be expanded as Taylor Series at $ x=0$ to obtain the sequence from the coefficients. However, when I try to solve it both using DSolve and NDSolve, I have no luck. With DSolve it just returns the request itself:

$ $ \text{DSolve}\left[\left\{0.5 x^5 y”(x)+(2. x+1) x^3 y'(x)+\left(1. x^3+x^2+x-1\right)y(x)+1=0,y(0)=1,y'(0)=1\right\},y,x\right]$ $

And with NDSolve I just receive errors and no equation:

Power::infy: Infinite expression 1/0.^5 encountered. Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0.`.

$ \text{NDSolve}\left[\left\{0.5 x^5 y”(x)+(2. x+1) x^3 y'(x)+\left(1. x^3+x^2+x-1\right)y(x)+1=0,y(0)=1,y'(0)=1\right\},y,\{x,0,1\}\right]$

How could I resolve this problem?

Solving a system of differential equations whose one of the coefficients is imported data

Suppose we have a coupled system of differential equations: \begin{equation} \frac{db}{dt}=(- \gamma_b -i\omega_b)b-i\frac{g}{2}p;\quad \frac{dp}{dt}=i\frac{g}{2}\Delta N(t) b-(\gamma_a+\gamma_b+2iJ)p. \end{equation} If $ \Delta N$ was fixed, the solution of the system would be like \begin{equation} \begin{pmatrix} b(t)\ p(t) \end{pmatrix}=\begin{pmatrix} a_{11}&a_{12}\ a_{21}&a_{22} \end{pmatrix}\begin{pmatrix} b(0)\ p(0) \end{pmatrix} \end{equation} Using the following code, I have found a $ 2\times 2$ matrix (called sol) whose entries are $ a_{ij}$ in the above equation:

rb=630;wb=75*10^6;g=0.63;ra=2.6*10^6;rm=3.6*10^6;J=6.3*10^7;DeltaN=0.164*10^5; m ={{-rb-I wb,-I g/2},{I g DeltaN/2,-(ra+rm+2 I J)}}; eigvec = Eigenvectors[m] // Transpose // Simplify; eigval = Eigenvalues[m] // Simplify; inv = Inverse[eigvec] // Simplify; v1 = eigval[[1]]; v2 = eigval[[2]]; sol = eigvec.{{E^(v1 t), 0}, {0, E^(v2 t) }}.inv;

If we suppose that $ p(0)=0$ , then one can easily plot $ |b(t)/b(0)|^2$ : simply plot $ a_{11}(t)$ . But the problem is that $ \Delta N$ is not fixed. It is a $ N\times 1$ matrix which I have obtained from another code written with Fortran and its type is data.txt. The elements of this file are calculated by assuming the time interval between each one is $ 0.001$ . That is, for $ t=0.001$ we have $ \Delta N_1$ , for $ t=0.002$ we have $ \Delta N_2$ , etc. But the time intervals are not included in the txt file.

One way that comes to my mind is this: Assuming we know the analytical form of solfor a fixed $ \Delta N$ , we set time, i.g., equal to $ 0.001$ and then substitute the first row of the txt file (I call it $ \Delta N_1$ ) into sol and find $ a_{11}$ . Then we raise time to $ 0.002$ , substitute $ \Delta N_1$ into sol, find $ a_{11}$ , and repeat the procedure to the last row of the txt file.

Now the question is this: how can I import the txt file to the code and do the procedure that I explained above to get some data like $ \{\{0.001,a11(0.001)\},\{0.002,a11(0.002)\},….\}$ where the first elements are time intervals and the second ones are $ a_{ij}$ corresponding to that particular time?

I had asked a similar question here enter link description here, but in that problem I did not have an external file with txt format.

I could not upload my txt file, so I write the first 10 elements if necessary:

0.164E+05

0.655E+05

0.146E+06

0.258E+06

0.400E+06

0.572E+06

0.776E+06

0.101E+07

0.129E+07

0.159E+07