Bottom Page

Thread Rating:
• 2 Vote(s) - 2 Average
• 1
• 2
• 3
• 4
• 5
 Fitting experimental data with differential equations : Optimization madoko Unladen Swallow Posts: 3 Threads: 1 Joined: Jan 2019 Reputation: 0 Likes received: 1 #1 Jan-14-2019, 03:52 PM Hello !!! I'm a physics student trying to solve an experimental problem in fluid dynamics and here is the issue I'm having. The data output of my experiment is a 2D trajectory ([X,Y] array). I also have a theoretical model in the form of 3 coupled differential equations, solved using Runge Kutta 4, which also gives me a 2D trajectory ([x,y] array). This model depends mainly on 3 constants (a,G,B) of unknown values. Here is an equivalent python script of my problem: ```import matplotlib.pyplot as plt import numpy as np XY = np.loadtxt("experimental_data.txt") plt.plot(XY[:,0],XY[:,0],label="experimental data") x,y = Solver_diff_eqt(a,B,G) #Solution of my model plt.plot(x,y,label="model") plt.legend() plt.show() ```Obviously my code is much more complicated, but this is the basic idea. So finally my question. I'm looking for a tool or an algorithm which could optimize the values of (a,G,B) by trying to minimize the difference between my XY graph and my xy graph. (I guess that a better way to say it is : "I've 2 distributions, one of which depends on 3 constants. I'm looking for a way to calculate these 3 constants so that the 2 distributions are the same.", but it's harder to understand) To give you some context on what I already tried. I tried doing it by hand with sliders but it's way too complicated. I tried to use scipy.optimize.curve_fit with an intelligent hack, but I'm obviously not intelligent enough . Now you know everything ... Can you please help me ?  Thanks You like this post micseydel Involuntary Spiderweb Collector Posts: 1,972 Threads: 49 Joined: Sep 2016 Reputation: 45 Likes received: 586 #2 Jan-14-2019, 11:37 PM I think your question is clear, but I'm not familiar with existing tools which do it. I think sympy is a good bet, though I'm not 100% sure you'll find what you're looking for there. Gribouillis  Posts: 1,656 Threads: 14 Joined: Jan 2018 Reputation: 154 Likes received: 391 #3 Jan-15-2019, 12:13 AM (This post was last modified: Jan-15-2019, 12:14 AM by Gribouillis. Edited 2 times in total.) You could define a functional that measures the distance between the two curves such as ```X, Y = ... def distance(a, B, G): x, y = Solver(a, B, G) return np.linalg.norm(y - Y) ```Then minimize this function of 3 variables with `scipy.optimize.minimize(distance, **args)`. micseydel likes this post madoko Unladen Swallow Posts: 3 Threads: 1 Joined: Jan 2019 Reputation: 0 Likes received: 1 #4 Jan-15-2019, 06:34 AM Thank you both for answering  (Jan-15-2019, 12:13 AM)Gribouillis Wrote: You could define a functional that measures the distance between the two curves such as ```X, Y = ... def distance(a, B, G): x, y = Solver(a, B, G) return np.linalg.norm(y - Y) ```Then minimize this function of 3 variables with `scipy.optimize.minimize(distance, **args)`. This looks promising Gribouillis ! It's going to take me some time implement because optimize.minimize has a lot of options and subtleties. I will post the solution to this problem if I'm successful. scidam  Posts: 514 Threads: 1 Joined: Mar 2018 Reputation: 69 Likes received: 69 #5 Jan-15-2019, 08:15 AM Hm..., this is quite interesting problem. Traditional approach to solving such problems consist in building iterative process converging to exact (or almost exact) solution (if this is possible). The approach proposed by Gribouillis is correct in general, but being applied in practice, it likely lead to very expensive and long computations. This is because Scipy's minimize function will likely try to compute gradient of the residual function `norm(y-Y)` that, in turn, depends on solution of the system of ODE (ordinary differential equations). Therefore, getting the gradient estimation will require a lot of computations. Another approach assumes the following steps: 1) Problem statement. Let we have (three ODE's as stated above) a system of ODEs and observations: Quote:dx/dt = F(x, y, p, a, B, G) dy/dt = G(x, y, p, a, B, G) dp/dt = H(x, y, p, a, B, G) O(t) = (x(t), y(t)).T # O - observation vector, .T means transposed,i.e. O is column vector Since a, B, G are constants the system above could be expanded as follows: Quote:dx/dt = F(x, y, p, a, B, G) dy/dt = G(x, y, p, a, B, G) dp/dt = H(x, y, p, a, B, G) da/dt = 0 dB/dt = 0 dG/dt = 0 O(t) = (x(t), y(t)).T This is a classical dynamical system, which input parameters, i.e. initial conditions of the ODE, are needed to be found. These parameters are: x0, y0, p0, a0, B0, G0. The last triple a0, B0, G0 doesn’t change in time (they are constants). So, we can rewrite the problem in vector-matrix notation as follows: Quote:dX/dt = Z(X) O(t) = M(X) + e(t) where X = (x, y, p, a, B, G).T; Z -- nonlinear function that returns a vector of size 6 and depends on X; O is a vector of available observations. It depends on time, and might be known only at specific time moments, e.g. t_k, k=1,...,N. In your case M is quite trivial, it returns just x and y. i.e. M = (X, X).T = (x, y).T; e(t) -- error of measurements vector (we don't need to know it, we will use it to state the optimization problem in least square sense). In general, Z(X) could be of non-linear and/or non-stationary (i.e. depends on time, Z = Z(X,t)) form. I hope that your ODEs is stationary one. So, if we have measurements (O(t_k)) on the trajectory, we need to find initial conditions X0=X(t0) that leads to the trajectory that, in turn, satisfies (in least square sense) to all of our observations. This means that the trajectory minimizes the sum(e(t_k)^2). Quote:sum e(t_k)^2 = sum (O(t_k) - M(X(X0, t_k))^2 => min by X0 This is non-linear (in general) optimization problem, but it could be solved in iterative way as a series of linear ones. Suppose that we have a good approximation of the solution, we denote it as R0. We can rewrite the system in variational form (full derivation of this I'm leaving out): Quote:dP/dt = A(R0) * P(t) dO_k = O(t_k) - O(R0, t_k) = H * P(t_k) + e(t_k) where P(t) = X(new_R0, t) - X(R0, t); X(new_R0, t) is a better estimation of the solution that we are trying to find; If we find P, we can get the better estimation from the current one, e.g. X(new_R0, t) = X(R0, t) + P or new_dR0 = R0 + P(t_0). O(t_k) -- our observations, we already know them; so, dO(t_k) is known vector of observation residuals; A(R0) is Jacobian matrix of Z(X), A(R0) = dZ/dX (when X=R0) (it doesn't depend on time if the system of ODEs is stationary); H is Jacobian matrix of the measurement function; in your case it would be very simple, H = [[1, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0]].T, since you're observing raw coordinates x and y, that are first two elements of X vector. The latter system is linear and it can be solved using matrix exponential. So, we have Quote:P(t_k) = EXP(A(R0)*(t_k-t0)) * P(t_0), and, finally, substituting the latter into equation of measurements we get an overdetermined system of linear equations: Quote:dO_k = H * EXP(A(R0)*(t_k-t0)) * P(t_0) + e(t_k), k = 1, ..., N (assume that we have N measurements) This is a system of linear algebraic equations (SLAE). It can be solved in least square sense, e.g. using Moore-Penrose transfrom/inverse (np.linalg.pinv); Once we get its solution P(t_0) we can improve our initial conditions: new_R0 = R0 + P(t_0) and repeat this process with new new_R0. Depending on choice of R0 this process could converge to the solution we are trying to find. This approach allows to get good estimation of initial conditions that leads to observed trajectory. Also, if you try to estimate not initial X(t0) but current X(t) state of the system, then you probably invent the Kalman filter... Hope that helps... You can read about linear/nonlinear dynamical systems, e.g. in classical work "Topics in mathematical system theory" by Kalman, Falb, Arbib. Gribouillis and snippsat like this post madoko Unladen Swallow Posts: 3 Threads: 1 Joined: Jan 2019 Reputation: 0 Likes received: 1 #6 Jan-16-2019, 12:26 AM First scidam, thank you so much for your answer. It's so complete and precise ! I will just add a few precisions: Yes my ODE system is indeed stationary. My initial conditions are known. I'm not sure I understand what you mean by: Quote:Suppose that we have a good approximation of the solution, we denote it as R0. The R0 here confuses me, especially in the case of this equation: Quote:dO_k = O(t_k) - O(R0, t_k) = H * P(t_k) + e(t_k) What you mean is that R0 = [a_approx,B_approx,G_approx], refining the approximation at each iteration ? Again thanks very much !  scidam  Posts: 514 Threads: 1 Joined: Mar 2018 Reputation: 69 Likes received: 69 #7 Jan-17-2019, 11:30 AM Because I didn't know that initial conditions are known I considered extended system, so R0 = [x0_approx, y0_appros, p0_approx, a_approx, B_approx, G_approx]; Therefore, dO_k -- residual between known measurement at time t_k and measurement at t_k obtained for trajectory started at (R0, t0). From now, we know that (x0, y0, p0) are exactly known (I wrote x0, y0, p0 because we have 3 coupled differential equations as you denoted that in the first post), so our iterative process should be slightly changed: There are several ways we can adopt our iterative process: 1) don't change x0, y0, p0 when iterating; apply changes to a, B, G only; 2) If you have measurement at t0, you can consider weighted least squares and assign very big weight to this measurement. Since measurements are raw coordinates big weight for t0-measurement will account that x0 and y0 are very important measurements and shouldn't be changed over iterations; Finally, everything is depending on the system of ODEs you are trying to solve. Is your system really hard to solve analytically? You are likely needing to construct problem-specific algorithm to get estimations for a, B, G. Have you some a priory information regarding a, B, G values...? If you have fast numerical procedure for solving the ODE, you can try, e.g. simulated annealing approach. At last, did you try the way proposed by Gribouillis? « Next Oldest | Next Newest »

Top Page

 Possibly Related Threads... Thread Author Replies Views Last Post Need Help solving second order differential equations SkewedZone 2 244 Jun-25-2019, 12:14 PM Last Post: SkewedZone Solve a system of non-linear equations in Python (scipy.optimize.fsolve) drudox 7 7,138 Aug-18-2018, 02:27 AM Last Post: scidam Python code optimization problem servanm 2 627 May-23-2018, 01:28 PM Last Post: servanm Fitting Lognormal Data Carolyn 3 3,535 May-21-2018, 12:10 AM Last Post: scidam I need a help to solve equations system alizo_as 1 539 May-04-2018, 04:51 PM Last Post: Gribouillis fsolve 3 equations with 3 unknowns. Can't get an answer. Aron313 6 1,074 Mar-18-2018, 04:35 PM Last Post: Aron313 Solving nth degree simultaneous equations hegdep 3 1,124 Sep-28-2017, 08:33 PM Last Post: nilamo

Forum Jump:

Users browsing this thread: 1 Guest(s)