6

The first step is to transform the second order equation to a set of two coupled first order equations. Define an auxiliary function $u(T) = \frac{dr(T)}{dT}$. This results in the system $$\begin{align} \frac{du}{dT} &= k-(1-\frac{5}{r})(3+\frac{2}{r^2}) \\ \frac{dr}{dT} &= u\\ \frac{d\phi}{dT} & = \frac{1}{r^2} \end{align} $$ Now you have a ...


5

I see at least one important problem. On the right hand side you have a term that looks like $$P_o \left( \frac{\dot{R}}{R} \right)^{3 \kappa}$$ This term is dimensionally inconsistent with the other terms in the brackets, which have dimensions of a pressure. This term should actually be $$P_o \left( \frac{R_o}{R} \right)^{3 \kappa}$$ The paper you link ...


5

You can convert this equation into a Poisson equation and use standard numerical methods to solve it like finite difference to obtain $\theta(x,y)$. If you take a divergence from both side: $$\nabla^{2} \theta = \frac{1}{k}\nabla \cdot \vec{v}$$ Since, you know the $\vec{v}$ you can approximate its divergence as: $$\nabla \cdot \vec{v} = \frac{\partial v_{...


4

You have to write your second order equation as a system of two first order equations. Let $y' = v$, then your equation $$ y'' + \omega^2 y = 0 $$ becomes $$ \begin{pmatrix} y' \\ v' \end{pmatrix} = \begin{pmatrix} v \\ -\omega^2 y \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -\omega^2 & 0 \end{pmatrix} \begin{pmatrix} y \\ v \end{pmatrix} $$ If you ...


4

You are starting from a uniform temperature and you have insulated boundary conditions; so there is no heat conduction occurring. Likewise, your initial mole fraction is also constant in $x$ so that the heat input is uniform along the length. So the fact that the temperature and mole fraction don't change as a function of $x$ is exactly what you should ...


4

Is there any out-of-the-shelf method I can use in this case to solve the equation without fully calculating the Jacobian? Using any integrator for stiff ODEs where the implicit equation is solved via Newton-Krylov methods will only require Jacobian-vector products. So something like DASKR, IDA, etc. can be made to do this.


3

You're approaching this the wrong way by using the product rule of differentiation. Rather, making use of the fact that $$ \frac{\partial u(x_i)}{\partial x} \approx \frac{u_{i+1/2}-u_{i-1/2}}{\Delta x}, $$ the first step in finding what you need to do is to recognize that $$ z(x_i)\frac{\partial u(x_i)}{\partial x} \approx z_i \frac{u_{i+1/2}-u_{i-1/2}...


3

It's hard to get around the allocations implicit to SymPy in this case. It wants to allocate the matrix, so the easiest thing to do would be, as you show, build individual scalar functions. But then composing those together can be a bit of a hassle, since you don't want to put them into an array since they are all different types and that would then ruin the ...


3

The parameter can be any type, so here I pass in a time-dependent function for p and use it in the differential equation: # Packages library(tidyverse) library(diffeqr) library(JuliaCall) diffeq_setup() # Drift function f <- function(u,p,t){ du1 = p(t) return(c(du1)) } # Diffusion function g <- function(u,p,t){ du1 = 0 # note that there is ...


2

The usual trick is to add more variables that represent the successive derivatives, as in the equations of motion pf physics written as a set of first order ODE of "2" variables: $$\begin{align}\dot{x}&=v\\ \dot{v}&=a(x,t)\end{align}$$ instead of a second order ODE of 1 variable $$\ddot{x}=a(x,t) \ .$$ The dot represents the time derivative. So, ...


2

In fact, your equation is a non-linear advection-diffusion. Due to the fact that your problem is time-dependent, it could be easily solved by finite-difference: $$\frac{\partial z}{\partial t} = -C (\sin(\omega t))^{m} x^{hm} n (\frac{\partial z}{\partial x})^{n-1} \frac{\partial^{2} z}{\partial x^{2}} -C (\sin(\omega t))^{m} x^{hm} h m x^{-1} (\frac{\...


2

Your cost function can also be written as $$ K = \int_0^{t_f} \left(\phi(t) - e^{-M^\top \tilde{D}\,M\,t} \hat{\phi}(0)\right)^\top \left(\phi(t) - e^{-M^\top \tilde{D}\,M\,t} \hat{\phi}(0)\right) dt. $$ When minimizing that cost function with respect to $\tilde{D}$ and $\hat{\phi}(0)$ it would be equivalent to minimizing the following cost function $$ K =...


2

I am a bit confused as to your characterization of constraints. Equation $(1)$ is not a constraint. It is the model that generated the time series data you are trying to fit. You then try to find the correct parameters $\tilde{D}$ that result in equation $(2)$ matching your time series as well as possible. I would formulate the problem as the following: ...


2

You're printing out the time steps. If you print out the solution, you should find the IC is initialized properly.


2

It seems like your goal is to get an accurate numerical solution for you differential equation, which likely does not require you to code your own ODE solver. In that case, it is likely more efficient for you to reframe your problem for use in an ODE solver in your programming environment of choice. A standard first step is to rewrite your equation as a ...


1

So the size of the problem doesn't necessarily stop you from using the full Jacobian if you store it in a compressed/sparse storage format, you also should be able to know which residuals are dependent on which unknowns which should make it cheaper as well. That said, just use Jacobian free Newton-Krylov (JFNK) solvers. FGMRES is a really strong linear ...


1

There's a really nice book by Li and Ito on the immersed interface method, which was designed to solve problems like yours. Chapter 2 describes 1D problems. Basically, the finite difference method works well when the coefficients are smooth, but when they're rapidly varying or have discontinuities things can go to hell very quickly. You can make it work by ...


1

First one should verify the numerical solution for $a(t)$ against the analytic solution. $ d_t a = 1/a + 1/a^2 = \frac{a+1}{a^2} $ Use $b=a+1$, then the ODE becomes $ \frac{(b-1)^2}{b} db = dt, $ from which we find the solution $ b^2/2 - 2 b - \ln(b) = t $ This relation (defined up to a constant to satisfy the initial conditions) is not invertible but ...


1

Using the given ODE, the boundary conditions $u(x)=0$ at $x=\pm1$, and the symmetry of the solution it is easy to write the exact analytic solution: $u(x)= a \sin(\pi x)$, for x $\in$ [1/4,1] $u(x)=- a \sin(\pi x)$, for x $\in$ [-1,-1/4] $u(x) = 1/\pi^2 + b \cos(\pi x)$, for x $\in$ [-1/4,1/4] The matching condition for the function derivative at x=1/4 ...


1

The integration curves have the overall form as shown in the figure below, they have vertical asymptotes at $y=0$, which certainly creates difficulties for a numerical solver integrating from $y=0$ (x does not matter). This does not necessarily mean the problem is ill posed; that would be the case if the integration curves diverge rapidly as the initial ...


1

Your problem is not well posed. Multiply both sides by $|y|$ to obtain the ODE in the form $$ |y(x)| y'(x) = A \sqrt{By(x)+C}. $$ If your initial condition is $y(0)=0$, then at the initial time you have $$ 0 y'(0) = A\sqrt{C}. $$ This means that if $A\sqrt{C}\neq 0$, no value of $y'(0)$ can satisfy the equation. In other words, your problem is not that ...


1

I have been thinking that it might be easier if one changes first the variables in the differential equation. That way one can bypass the function $h(t)$ and deal with fewer functions. Since $$\hat{y}(t) = \hat{y}_{\theta}(t) = \hat{y}(t \,| \, \theta) = \theta \, \sqrt{2g} \, \sqrt{h(t)}$$ change the dependent variable $$\hat{y} = \theta \, \sqrt{2g} \, \...


1

I don't see how two equations give $z$ as output. Nevertheless, your sequence of computations looks reasonable, except I would combine steps one and two into: Simultaneously solve for $x$ and the sensitivities $z$. This is a extended ODE system that you could throw at a built-in MATLAB ODE solver: $$ \begin{bmatrix} x'(t,\hat{\theta}) \\ z'(t,\hat{\theta}...


1

Weighted least squares is a method to induce a-priori knowledge into the ordinary least-squares procedure. The basic example is a simple scaling of the dimensions, such that each predictor variable has about the same variance (which is up to a factor also what you found in your reference). In your case, you would use that if the parameters $y$, $\dot y$, ...


1

This is a system of first order differential equations, not second order. It models the geodesics in Schwarzchield geometry. In other words, this system represents the general relativistic motion of a test particle in static spherically symmetric gravitational field. In general, there is a third equation for how coordinate time is related to proper time. ...


1

In this code, you are trying to solve this 1D heat transfer equation: $$\frac{\partial u}{\partial t} = \alpha \frac{\partial^{2} u}{\partial x^{2}}$$ If you use discrete Fourier transform (DFT) on $u(x,t)$ and discretized $x$ values of $x_{i}$ where $0 \leq i \leq N$ and $x \in [0,L]$ and $x_{i} = \frac{i}{N}L$ and $k_{i} = \frac{2 \pi i}{L}$, you have: ...


1

I think the problem with initial conditions is more related with the sum of very small numbers with other that are only small (I suspect the problem resides in terms like gam^8, since gam is small). The resulting numeric inaccuracy will bring problems. I've played around with other initial conditions and integration times and in octave (that is similar but ...


1

Take the first expression and start to reduce the $x$ values to $x_i$, \begin{align} \frac{u_{i+1}^n - u_i^n}{x_{i+\frac{1}{2}}} - \frac{u_i^n - u_{i-1}^n}{x_{i-\frac{1}{2}}} &= \frac{(x_i-\frac12Δx)(u_{i+1}^n - u_i^n) - (x_i+\frac12Δx)(u_i^n - u_{i-1}^n)}{x_{i-\frac{1}{2}}x_{i+\frac{1}{2}}} \\ &=\frac{x_i}{x_i^2-\frac14Δx^2}(u_{i+1}^n - 2u_i^n + u_{...


Only top voted, non community-wiki answers of a minimum length are eligible