2,189 questions
0
votes
0
answers
31
views
Solving differential equations with ode solvers and finite difference in Matlab
I have implemented following equations both with the Matlab own ODE solvers and with a very simple finite difference scheme. The latter works properly, while the ODE code does not produce suitable ...
2
votes
1
answer
137
views
Problem with mismatched length when using a mask
I'm writing a code and I have a function that calculates the values that are not fulfilling a condition with the values that are fulfilling the condition, but I'm having a lot of trouble with managing ...
0
votes
1
answer
63
views
Numerical instability in forward-backward algorithm for Hidden Markov Models
I am implementing the forward algorithm for Hidden Markov Models (see below for the algorithm). To prevent over/underflow, I work with log-probabilities instead, and use the log-sum-exp trick to ...
0
votes
0
answers
16
views
How to accurately draw a Bode diagram of a control system in MATLAB in the case of its matrix being ill-conditioned?
I'm trying to draw a bodeplot of a discrete contorl system with the "ss" and "bodeplot" function:
% X(k+1) = AX(k)+Bu(k)
% y(k) = CX(k)+Du(k)
sys_d = ss(m_A, m_B, m_C, m_D, Ts);
h1 ...
3
votes
2
answers
118
views
Using SciPy minimize to solve the Catenary problem
I am attempting to use SciPy optimization to find the function that minimizes the potential energy of a hanging rope between 2 points, the catenary problem. The objective function is defined by
The ...
2
votes
2
answers
137
views
Algorithm for calculating fraction from decimal with limits
The algorithm for calculating a fraction from a decimal number is known and there are enough examples for that. In my case I need to find a fraction for a given decimal where A = N/D is optimized but ...
2
votes
0
answers
76
views
Adams-Bashforth coefficients
I implemented a first draft of the Adams-Bashforth multistep-solver in C++. To make this method work, you need coefficients. In Wikipedia these are the bj in this formula:
y_n+1 = y_n + h * sum_over_j(...
1
vote
2
answers
151
views
Problem with least squares rational approximation to `asin(x)+ sqrt(1-x^2)` in [3,1] form
I'm trying to generate a decent [3,1] rational least squares polynomial approximation to asin(x)+sqrt(1-x^2) on [0,1] and failing dismally :(
The problem is that it has a pole for this particular ...
1
vote
0
answers
43
views
Solving Nonlinear Boundary Value Problem with Constraints
I want to solve the following boundary value problem
where S, k_i, x_f, α_1, and θ are known parameters. We are trying to solve for h(x), p, and θ_d.
My idea was to use finite differences to create a ...
0
votes
0
answers
25
views
Numerical solve in NetLogo
We are using a mass balance model (Dynamic Energy Budget Theory based), and know that for 2 (maybe 3) parameters, we will need to use numerical solving methods to calculate fluxes (it says this in the ...
1
vote
0
answers
150
views
A more efficient and robust algorithm for minimizing an objective with many local minima
Background
Background
I've developed a Python script that calculates the time taken for a ray of light to propagate between two points in a complex medium. Now, I want to extend this script to solve ...
2
votes
0
answers
56
views
Locating a real-valued function in a 2-dimensional matrix [closed]
I need help solving the problem of locating a chosen function in a 2-dimensional matrix. Let me explain what I mean precisely. I have two 1-dimensional real grids (arrays) of x(i) and y(j), where i, j ...
0
votes
0
answers
31
views
My code for LU factorization in matlab is not working
I have typed this code to perform the LU factorization for A=xb when we are given matrix A and $x=(1,1,1)^T$
. My code is below
function q1()
% Define x_ex
x_ex = [1; 1; 1];
% Define ...
0
votes
0
answers
29
views
Integrating a 2nd ODE that contain a first derivative squared term
The 2nd ODE that I want to solve numerically is:
xddot = ( 0.5 * rho * A * Cl * ( xdot^2 + vy^2) ) / m
ODE
The vy term is a sinusoidal function and does not need to be solved numerically.
I have ...
0
votes
0
answers
30
views
Truncated Newton method with preconditioning on MATLAB
I'm trying to implement TN with preconditioning on MATLAB, and evaluate it for different functions, with different starting points. I'm also trying to confront it with TN without preconditioning.
What ...
3
votes
1
answer
83
views
MSVC FP library float erfcf(x) unexpectedly slow - 2x slower than double erfc(x). Why?
Can someone please verify that the slowness that I see with the MS library implementation of erfcf(x) is reproducible elsewhere. It might just be a benchmarking quirk. Also any suggestions for command ...
0
votes
0
answers
78
views
Finding Local/Global Maxima of high-dimensional function in Mathematica
I am trying to find the local/global maxima of a function. This function is computed by taking the trace of a 64x64 matrix and summing over a bunch of indices. With my current implementation I have 90 ...
0
votes
0
answers
37
views
The following code is to solve a hydrodynamic stability problem. It is a shooting method implemented in Octave software
I have the following code written in Octave. It is supposed to solve a coupled differential equation (DE) dealing with hydrodynamic stability analysis. The code works fine without imaginary numbers ...
0
votes
1
answer
69
views
Plotting a nonlinear function in R [closed]
I am currently struggling to plot the next function in R.
(F1 / F2^m * (R1 - F1)^((1 - s) / s)) - (m * (((R1 - F1))^(1 / s) + ((R2 - F2))^(1 / s)) / (F1^m + F2^m)) = 0
Where I want to plot F2 as a ...
0
votes
0
answers
30
views
Calibration, Sensitivity Analysis, Numerical dispersion And Model Verification In An Ecological Model With 15 state variables
I am a student assigned to a project on improving an ecological model. This model is 1-D model consisting of a dispersive-advective component and pelagic-reaction component which was developed in R. ...
0
votes
0
answers
40
views
Determine divergence points in root search of a non-continuous function
I need to implement a root-searching algorithm for a function, that can possibly have some (non-removable) divergence points (DPs)/singularities, e.g. f(x) = x**2 * np.tan(x).
The problem is, that the ...
0
votes
0
answers
53
views
How to optimize the function below so that it might run under 1 second
Below is a code that makes a network where each node of the network is a Stuart-Landau oscillator. The equations of evolution are coupled and also depend on it's neighbors. The code runs 10^6 ...
0
votes
1
answer
41
views
SciPy's solve_bvp going solving for the wrong parameter values
I am solving a ordinary differential equation with scipy's solve_bvp to get a family of eigenstates phi_n and eigenvalues omega_n. The equation looks something like
y'' = -(omega_n - f(x)) ** 2 * y,
...
2
votes
0
answers
39
views
cupy.linalg.solve for positive definite matrix?
It seems like cupy.linalg.solve doesn't have an option for me to solve linear system Ax=b assumingA is positive definite?
I am looking for something like scipy.linalg.solve where one can actually tell ...
0
votes
1
answer
147
views
Using `solve_ivp` and `LSODA` to solve a complex ODE
I am trying to solve a complex system of differential equations. The equations are stiff so I need to need to use a method which can handle both complex ODEs and stiffness switching. I have landed on ...
1
vote
0
answers
105
views
TOMS 748 - "the bounds given must bracket a single root"?
I have a cubic equation f(x)=0 which I need to solve numerically to within some tolerance. I only need one of the roots.
I already have values a and b such that f(a) and f(b) have opposite signs; i.e. ...
0
votes
3
answers
73
views
python and problems with passing array to a function
I need to solve some ordinary differential equation dy/dx = f(x) = x^2 ln(x) and to proceed I create the array xpt between the limits 0. <= xpt <= 2. Because I have to be careful at xpt = 0, I ...
3
votes
1
answer
122
views
Inaccurate Eigenvectors from LAPACK's ZHEEVR routine
I was trying to replace Scipy's eigh routine in my cython code with C-level LAPACK. Essentially I want to compute eigenvectors corresponding to the largest and smallest eigenvalues of a hermitian ...
0
votes
0
answers
54
views
Generate Hopf bifurcation diagram for 3D system in Matlab
Hello I want to plot the Hopf bifurcation diagram for the 3D system.
I have this code but I'm getting nothing other than straight lines.
How would I fix this? Any help is appreciated. The parameter ...
2
votes
1
answer
97
views
Optimize this Python code that involves matrix inversion
I have this line of code that involves a matrix inversion:
X = A @ B @ np.linalg.pinv(S)
A is an n by n matrix, B is an n by m matrix, and S is an m by m matrix. m is smaller than n but usually not ...
3
votes
2
answers
189
views
How to solve complex equations numerically in Python?
Originally, I had the following equation:
2mgv×sin(\alpha) = CdA×\rho(v^2 + v_{wind}^2 + 2vv_{wind}cos(\phi))^(3/2)
which I could express as the following non-linear equation:
(K × v)^(2/3) = v^2 + ...
1
vote
1
answer
69
views
CFD Navier-Stokes Equation in Julia
so below is the code snippet i have for Julia:
using Printf
using Plots
gr()
# Parameters definition
Ny = 21 # Number of points in y-direction
Aspect_ratio = 10 # Ratio of ...
1
vote
0
answers
36
views
RK4 acoustic wave equation
I am trying to modify a code to include a runge kutta time stepping scheme. The original code was not written by me but my motive is to eliminate the pygame library and make it as easy as possible ...
2
votes
2
answers
101
views
Creating a short, d-dense list of points on the sphere S^4 in Python
I am a mathematician and this is my first time using Stack Overflow, sorry if the question is not adequate or there is some better place to ask this. I would like to know if there is a standard way ...
1
vote
1
answer
96
views
What is going on behind this numerical integration?
I am calculating the value of the integral
∫₀¹ ∫₁² (r - u)r drdu
I would like to understand how the following error is handled in R.
Looking around the forum, I realized that a vectorization of the ...
0
votes
0
answers
52
views
How to do a matrix optimization by matlab?
I want to solve this optimization problem (which is a simplified demo) by matlab. The optimization variables are A0,B0 which are 4 × 4 symmetric matrices. The constraint is -I≤ A0,B0≤ I. The objective ...
2
votes
1
answer
107
views
Problems with a simple numerical integration in Python
I am trying to solve a simple integral in Python using different methods for numerical integration. However, when I try to calculate the integrals, I obtain wrong results for small values of my x ...
0
votes
1
answer
88
views
Runge-Kutta Behaving Strangely When Initial Value Is Negative
I have the following RK4 algorithm in C++ for solving first-order differential equations:
using ODE_Function = std::function<double/*dy/dx*/(double/*x*/, double/*y*/)>;
template<int length&...
0
votes
0
answers
54
views
2D-FDM of heat exchanger
I am an undergraduate student majoring in Mechanical Engineering.
I am currently attempting to compare the e-NTU method and the Energy Balance Equation for a 2D-Laminar counter flow heat exchanger. ...
0
votes
0
answers
35
views
Krylov update after column change
(cross-post from https://math.stackexchange.com/questions/4928836/krylov-update-after-column-change)
Having a Krylov space of A
[V, H] = Arnoldi(A, v1, m),
how to update V and H after adding vector ...
0
votes
0
answers
34
views
calculate wave speed using Crank-Nicolson Finite diffrence sheme
im studying travelling waves for multidimension reaction diffusion equation so i start with 2D Fisher KPP to do , i use implicit Crank-Nicolson Finite diffrence sheme
here the code:
r = 1.0;
D = 1.0;...
0
votes
1
answer
228
views
How do I use numpy.gradient() with coordinate arrays?
I have three 1D arrays representing some spatial data (u) over a rectangular region of space: the first array contains the x coordinates, the second contains the y coordinates, and the third contains ...
0
votes
1
answer
78
views
Verify Solution Satisfies a PDE: FiPy
Say I have a solution to the following implementation of the Heat Equation using FiPy:
from fipy import CellVariable, FaceVariable, TransientTerm, DiffusionTerm, Grid1D
from fipy.tools import numerix
...
0
votes
0
answers
35
views
ADI method in matlab
this is code of 2D reaction diffusion equation using ADI method the code works well when i use the exponential
this is the code :
r = 1;
D = 1;
h = 0.0025;
Sx = 40.0;
Lx = -20.0;
Sy = 0.5;
Ly = -0.5;
...
1
vote
1
answer
59
views
Learning OOP in MATLAB - what kind of method/function/routine to use to construct a problem-specific vector?
I am working on a project in MATLAB that is essentially a linear system solver for a physical problem with lots of parameters and this is my first time trying OOP. Before, I had lots of MATLAB scripts ...
0
votes
0
answers
20
views
wave spped of propagation of solution by Matlab
ipm exploroinf DFisher KPP equation and diffrent numerical methods to study its behaviour specifically travelling waves and speed of propagation
here is the code in Matlab where iuse finite ...
0
votes
0
answers
24
views
Slope graph and solution
i am doing finite difference scheme to solve Fisher kpp equation :$u_{t}=u_{xx}+u_{yy}+u-u^2 and with initial condition u(x,y,0)=1 for (x,y)=(0,0)and 0 other wise and u=0 at boundary ,also plot ...
0
votes
0
answers
13
views
finite difference scheme and visualisation of the solution
i try to solver $u_{t}=u_{xx}+u_{yy}+u-u^2 in(-a,a)\times(-b,b)\times(0,T)with initial condition u_0=1 if(x,y)=0and 0 other wise and u=0 in$\partial\omega$ where $\omega=[-a;a]\times[-b;b], i do that ...
2
votes
0
answers
98
views
C++ Hermite interpolation - Generate coefficients & value
I'm trying to fix my Hermite interpolation program, which is supposed to output the coefficients and value at a point, but it seems the results are wrong. It currently looks like this:
#include <...
0
votes
1
answer
379
views
How to quickly measure the rise time of a signal represented inside a numpy array
I have a numpy array called voltage which is supposed to go from 0 to 1, in 1ps time steps. Thing is some times it may not go from 0 to 1, it might rise from zero and come back to 1.
I want to measure ...