Using the RungeKutta integration method in a system
h=0.005;
x = 0:h:40;
y = zeros(1,length(x));
y(1) = 0;
F_xy = ;
for i=1:(length(x)1)
k_1 = F_xy(x(i),y(i));
k_2 = F_xy(x(i)+0.5*h,y(i)+0.5*h*k_1);
k_3 = F_xy((x(i)+0.5*h),(y(i)+0.5*h*k_2));
k_4 = F_xy((x(i)+h),(y(i)+k_3*h));
y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;
end
I have the following code, I think it's right. I know there's parts missing on the F_xy because this is my follow up question.
I have dx/dt = = −x(2 − y) with t_0 = 0, x(t_0) = 1
and dy/dt = y(1 − 2x) with t_0 = 0, y(t_0) = 2.
My question is that I don't know how to get these equations in to the code. All help appreciated
1 answer

Is
F_xy
your derivative function?If so, simply write it as a helper function or function handle. For example,
F_xy=@(x,y)[x*(2y);y*(12*x)];
Also note that your
k_1, k_2, k_3, k_4, y(i)
are all twodimensional. You need to resize youry
and rewrite the indices in your iterating steps accordingly.
See also questions close to this topic

Simulink error: Initial conditions solve failed to converge.
I have the circuit in this figure:
I get the following error: Initial conditions solve failed to converge. Nonlinear solver: Linear Algebra error. Failed to solve using iteration matrix. The model may not give enough information to make it possible to solve for values of some of its variables. Specific advice is given below.
all components and nodal across variables involved Tie variable x to a definite value, for example by connecting an appropriate domain reference block.
I have looked at the other questions on this on MATLAB answers (did not see any on here) and these solutions did not work:
 Changing the SPS block to filter input.
 One question was missing ground but I have that.
Any idea what is wrong?

plot multiple Occupancy grid in one plot in matlab
How can I put multiple Occupancy Grid in one plot in matlab?
I tried to put them like an image using subplot but it is not working.

MATLAB 3D surfaces
in MATLAB, how do u plot a 3D surface using fsurf. I need to plot ((exp(t))sint,(exp(t))cost,1) but the plot will not load. Do I have the vector in the wrong format?

Interpolation of scattered data: scipy.interpolate.griddata appears to have prohibitive overhead
I'm currently porting some code from Matlab to Python. Basically, it's a numeric PDE solver in a 2D domain, which solves iteratively in slices of constant x. The input I'm given is scattered in 2D, which must be interpolated onto a regularlyspaced 1D grid at each slice. Here is the snippet that invokes the interpolator, in Matlab:
v_interpolator = scatteredInterpolant(x_scatt, y_scatt, v_scatt, 'linear', 'none'); yi = y_grid; % y_grid has about a thousand elements for i = 1:length(x_grid) xi = x_grid(i)*ones(length(y_grid)) v = v_interpolator(xi, yi); ...
Here it is in Python. Python does not have the equivalent of Matlab's
scatteredInterpolant
function. I wrote this class to make the code look as similar as possible to Matlab's builtin:from scipy.interpolate import griddata class ScatteredInterpolant: def __init__(self, x, y, z): obj.x = x obj.y = y obj.v = v def __call__(xi, yi): return griddata((obj.x, obj.y), obj.v, (xi[None,:],yi[:,None]), method='linear') ... v_interpolator = ScatteredInterpolant(x_scatt, y_scatt, v_scatt) yi = y_grid # y_grid has about a thousand elements for i in range(len(x_grid)): xi = x_grid[i] v = v_interpolator(xi, yi) ...
After timing these two different versions, I find that the interpolation is prohibitively slow in Python, for
y_grid
of size 1000: Matlab: each interpolation takes ~3 ms
 Python: each interpolation takes ~1.6 s, more than 500 times worse!
At some point it occurred to me to do all the interpolation in one step rather than inside the forloop, to create a single 2D grid of interpolated values. Then, inside the loop, I could just index the slice that I need.
# do all interpolation _once_, incurring the 1.6 s overhead only once v_grid = v_interpolator(x_grid, y_grid) for i in range(len(x_grid)) # index into grid v = v_grid[i, :] ...
For my purposes, this is impractical, as the entire grid will not fit in memory. Nevertheless, I tested this approach, and discovered that there is actually an overhead to calling
griddata
: as I increase the number of points in the xdirection of the grid that I interpolate onto, the time doesn't really change until the grid size becomes appreciably large. This overhead is clearly seen as the red line in the graph below: it's a minimum cost under which the timing will never go below, no matter how small the interpolation grid is.This doesn't happen in Matlab's builtin
scatteredInterpolant
, which readily scales with the number of grid points. My question is, why is Matlab's implementation so much faster, and/or how can I achieve similar performance in Python?From this post, it appears that I will need to dig more into what's going on behind the scenes in
griddata
. But I'm also curious as to why this isn't necessary in Matlab. 
Programming interpolation spline with R
I am stuck at programming interpolation spline. I did some math, which I have to code up:
Fi(x) = A1 + (B1 / x+C1) => Fi1(x1)=f1, Fi1(x0)=f0, F1`(x0) = g0
Then
A1 + (B1 / x1+C1) = f1 A1 + (B1 / x0+C1) = f0 B / (x0+C2)^2 = g0
Next I I need A1, B1, C1 (calculate somehow) Then calculate:
g1 = B1/(x1+c1)^2 = Fi1`(x1) = Fi2`(x1)
Please anybody has any knowledge how to start and calculate anything?
This is first deduction I got

Why Tensorflow can not solve this kinds of problem?
data.py is used to generation sample data. test.py is main code for parameters solving.
Objects movement in 3D space can be described by 4x4 Matrix, and this kinds of Matrix can be represented by 6 parameters. I have a group of Objects position data, and another group of final position data after Objects transformation. I want to use tensorflow for numeric computation of these 6 parameters. It shows a error of 'model trainning: None' when I run following code(test.py). Anybody can help me to debug this? Is it a limitation of tensorflow for this kinds matrix computation problem?
data.py
# coding: utf8 # In[1]: import numpy as np import tensorflow import pandas as pd # In[2]: def sample_points(low=0,high=100,N=10000): '''how to stack 1darray to 2darray''' x = np.random.randint(low,high,N) y = x*2 z = np.random.randint(low,high,N) L = np.full((1,N),1,dtype=np.int) return np.array([x,y,z,[1 for i in range(N)]]) # In[3]: xyzL_h=sample_points() print(xyzL_h) # In[4]: def save_points(xyzL,filename='xyz.txt',homogenous=False): ''' save xyz to ?.csv file, or xyz1 to ?.csv file ''' xyzL_v = np.transpose(xyzL,(1,0)) #print(xyzL_v.shape) df = pd.DataFrame(xyzL_v) df.columns=['x','y','z','1'] if not homogenous: df.iloc[:,0:3].to_csv(filename,index=False,header=False) else: df.to_csv(filename,index=False,header=False) # In[5]: save_points(xyzL_h,filename='raw_xyz.txt') # In[6]: import matplotlib.pyplot as plt get_ipython().magic('matplotlib inline') # In[7]: def show_data(xyzL_h): plt.plot(xyzL_h[0],xyzL_h[1],'r') plt.plot(xyzL_h[0],xyzL_h[2],'b+') # In[8]: show_data(xyzL_h) # In[9]: class Screw(): ''' based on screw thoery of screw motion and transformation matrix ''' def __init__(self,lib): self._lib = lib self._p = None self._v = None self._w = None self._seta = 0 self._m = None def dot(self,x,y): if self._lib.__name__ is 'numpy': return self._lib.dot(x,y) if self._lib.__name__ is 'tensorflow': #print("rank of y: ",y.shape,len(y.shape)) if len(y.shape)==1: y = tensorflow.reshape(y,(1,1)) return self._lib.matmul(x,y) def cross(self,x,y): ''' cross product of two vertors ''' x_n = self._lib.reshape(x,(1,3)) y_n = self._lib.reshape(y,(1,3)) return self._lib.cross(x_n,y_n) def transform(self,cloud): ''' cloud with shape of [4,None] ''' return self.dot(self._m,cloud) def I(self,n): ''' unit matrix ''' if self._lib.__name__ == 'numpy': return self._lib.identity(n) elif self._lib.__name__ == 'tensorflow': return self._lib.matrix_diag([1.0 for i in range(n)]) def w_from_rad(self,angle): ''' define vector by sphere coordinate first of angle(alpha): is radian angle between direction vector and Z axis second of angle(beta): is radian angle between XOY projection of vector and X axis ''' lib = self._lib alpha,beta = angle[0],angle[1] k = lib.cos(alpha) i = lib.sin(alpha)*lib.cos(beta) j = lib.sin(alpha)*lib.sin(beta) if lib.__name__ == 'numpy': ijk = lib.array([i,j,k]) elif lib.__name__ == 'tensorflow': ijk = lib.stack([i,j,k]) self._w = ijk return ijk def wx(self,w): ''' #skew symmetric matrix of twist component[0:3] w ''' if self._lib.__name__ is 'numpy': m = self._lib.zeros((3,3)) m[0,1]=w[2] m[0,2]=w[1] m[1,0]=w[2] m[1,2]=w[0] m[2,0]=w[1] m[2,1]=w[0] elif self._lib.__name__ is 'tensorflow': m = self._lib.stack([[0,w[2],w[1]],[w[2],0,w[0]],[w[1],w[0],0]]) # m = self._lib.reshape(m,(3,3)) return m def v(self,w,p): ''' # twist component[3:6] v ''' _v = self.dot(self.wx(w),p) self._v = _v return _v def tm(self,w,p,seta): ''' transformation matrix out of screw theory this matrix is used for right multiplication operation of point cloud ''' lib = self._lib if lib.__name__ is 'numpy': w = lib.array(w) p = lib.array(p) if lib.__name__ is 'tensorflow': w = lib.to_float(lib.Variable(w)) p = lib.to_float(lib.Variable(p)) _v = self.v(w,p) _wx = self.wx(w) self._w = w self._v = _v self._p = p m33 = self.I(3)+_wx*lib.sin(seta)+self.dot(_wx,_wx)*(1lib.cos(seta)) mm = self.I(3)m33 #+w*self.dot(w,_v)*seta #'w*np.dot(w,_v)*seta' should be 0, when no translation movement vv = self._lib.reshape(self.cross(w,_v),(1,1)) t = self.dot(mm,vv)[:,0] if lib.__name__ is 'numpy': m44 = self.I(4) m44[0:3,0:3]=m33 m44[0:3,3]=t self._m = m44#lib.transpose(m44,(1,0)) elif lib.__name__ is 'tensorflow': #print('t:',t) m34 = self._lib.stack([m33[:,0],m33[:,1],m33[:,2],t],axis=1) m44 = self._lib.stack([m34[0],m34[1],m34[2],[0.0,0.0,0.0,1.0]]) self._m = m44#lib.transpose(m44,perm=[1,0]) return self._m # In[10]: screw = Screw(np) q,seta=[0.0,100.0,0.0],np.pi/2 w =screw.w_from_rad([np.pi/4,np.pi/4]) res = screw.tm(w,q,seta) # In[11]: print(res) # In[12]: class CloudProcess(): def __init__(self,cloud,screw): self._cloud = cloud self._matrix = screw def transform(self): new_cloud = self._matrix.transform(self._cloud) return new_cloud # In[13]: op = CloudProcess(xyzL_h,screw) xyzL_h_new = op.transform() save_points(xyzL_h_new,filename='new_xyz.txt') # In[14]: screw = Screw(tensorflow) clouds = tensorflow.placeholder(dtype=tensorflow.float32,shape=[4,None],name='XYZL') p = tensorflow.Variable([0,100,0],name='PointOfAxis',dtype=tensorflow.float32) angle = tensorflow.Variable([np.pi/4,np.pi/4],name='DirectionAngle',dtype=tensorflow.float32) seta = tensorflow.constant(np.pi/2,dtype=tensorflow.float32) w = screw.w_from_rad(angle) m = screw.tm(w,p,seta) op = CloudProcess(clouds,screw) xyzL_new = op.transform() # In[15]: init = tensorflow.global_variables_initializer() with tensorflow.Session() as sess: sess.run(init) res=sess.run(xyzL_new,feed_dict={clouds:xyzL_h}) save_points(res,filename='new_xyz_tf.txt') # In[ ]:
test.py
# coding: utf8 # In[1]: import tensorflow import pandas as pd import numpy as np # In[2]: def load_pointXYZ(filename): df = pd.read_csv(filename,header=None) df['1']=1 df.columns = ['X','Y','Z','1'] return df # In[3]: df1 = load_pointXYZ('raw_xyz.txt') df2 = load_pointXYZ('new_xyz_tf.txt') # In[4]: X = df1.as_matrix().T Y = df2.as_matrix().T print(X) print(Y) # In[5]: class Screw(): ''' based on screw thoery of screw motion and transformation matrix ''' def __init__(self,lib): self._lib = lib self._p = None self._v = None self._w = None self._seta = 0 self._m = None def dot(self,x,y): if self._lib.__name__ is 'numpy': return self._lib.dot(x,y) if self._lib.__name__ is 'tensorflow': #print("rank of y: ",y.shape,len(y.shape)) if len(y.shape)==1: y = tensorflow.reshape(y,(1,1)) return self._lib.matmul(x,y) def cross(self,x,y): ''' cross product of two vertors ''' x_n = self._lib.reshape(x,(1,3)) y_n = self._lib.reshape(y,(1,3)) return self._lib.cross(x_n,y_n) def transform(self,cloud): ''' cloud with shape of [4,None] ''' return self.dot(self._m,cloud) def I(self,n): ''' unit matrix ''' if self._lib.__name__ == 'numpy': return self._lib.identity(n) elif self._lib.__name__ == 'tensorflow': return self._lib.matrix_diag([1.0 for i in range(n)]) def w_from_rad(self,angle): ''' define vector by sphere coordinate first of angle(alpha): is radian angle between direction vector and Z axis second of angle(beta): is radian angle between XOY projection of vector and X axis ''' lib = self._lib alpha,beta = angle[0],angle[1] k = lib.cos(alpha) i = lib.sin(alpha)*lib.cos(beta) j = lib.sin(alpha)*lib.sin(beta) if lib.__name__ == 'numpy': ijk = lib.array([i,j,k]) elif lib.__name__ == 'tensorflow': ijk = lib.stack([i,j,k]) self._w = ijk return ijk def wx(self,w): ''' #skew symmetric matrix of twist component[0:3] w ''' if self._lib.__name__ is 'numpy': m = self._lib.zeros((3,3)) m[0,1]=w[2] m[0,2]=w[1] m[1,0]=w[2] m[1,2]=w[0] m[2,0]=w[1] m[2,1]=w[0] elif self._lib.__name__ is 'tensorflow': m = self._lib.stack([[0,w[2],w[1]],[w[2],0,w[0]],[w[1],w[0],0]]) # m = self._lib.reshape(m,(3,3)) return m def v(self,w,p): ''' # twist component[3:6] v ''' _v = self.dot(self.wx(w),p) self._v = _v return _v def tm(self,w,p,seta): ''' transformation matrix out of screw theory this matrix is used for right multiplication operation of point cloud ''' lib = self._lib if lib.__name__ is 'numpy': w = lib.array(w) p = lib.array(p) if lib.__name__ is 'tensorflow': w = lib.to_float(lib.Variable(w)) p = lib.to_float(lib.Variable(p)) _v = self.v(w,p) _wx = self.wx(w) self._w = w self._v = _v self._p = p m33 = self.I(3)+_wx*lib.sin(seta)+self.dot(_wx,_wx)*(1lib.cos(seta)) mm = self.I(3)m33 #+w*self.dot(w,_v)*seta #'w*np.dot(w,_v)*seta' should be 0, when no translation movement vv = self._lib.reshape(self.cross(w,_v),(1,1)) t = self.dot(mm,vv)[:,0] if lib.__name__ is 'numpy': m44 = self.I(4) m44[0:3,0:3]=m33 m44[0:3,3]=t self._m = m44#lib.transpose(m44,(1,0)) elif lib.__name__ is 'tensorflow': m34 = self._lib.stack([m33[:,0],m33[:,1],m33[:,2],t],axis=1) m44 = self._lib.stack([m34[0],m34[1],m34[2],[0.0,0.0,0.0,1.0]]) self._m = m44#lib.transpose(m44,perm=[1,0]) return self._m # In[6]: class CloudProcess(): def __init__(self,cloud,screw): self._cloud = cloud self._matrix = screw def transform(self): new_cloud = self._matrix.transform(self._cloud) return new_cloud # In[7]: screw = Screw(tensorflow) cld_1 = tensorflow.placeholder(dtype=tensorflow.float32,shape=[4,None],name='XYZL') cld_2 = tensorflow.placeholder(dtype=tensorflow.float32,shape=[4,None],name='XYZL') p = tensorflow.Variable([0,50,2],name='PointOfAxis',dtype=tensorflow.float32) angle = tensorflow.Variable([0.3,0.2],name='DirectionAngle',dtype=tensorflow.float32) seta = tensorflow.constant(np.pi/2,dtype=tensorflow.float32) w = screw.w_from_rad(angle) m = screw.tm(w,p,seta) op = CloudProcess(cld_1,screw) xyzL_new = op.transform() loss = tensorflow.reduce_sum(tensorflow.abs(xyzL_newcld_2)) optimizer = tensorflow.train.GradientDescentOptimizer(learning_rate=0.1) model = optimizer.minimize(loss) init = tensorflow.global_variables_initializer() with tensorflow.Session() as sess: sess.run(init) for i in range(10): res = sess.run(model,feed_dict={cld_1:X,cld_2:Y}) print('model tranning:',res) res = sess.run(loss,feed_dict={cld_1:X,cld_2:Y}) print('loss:',res) print(angle.eval()) print(p.eval() )

Does anybody know why integrate in R breaks for large numbers in a beta?
I have this joint beta distribution I'm trying to integrate over.
alpha1 = 400 beta1 = 26000 alpha2 = 410 beta2 = 26000 integrate(function(y) { sapply(y, function(y) {integrate(function(x) dbeta(x, shape1 = alpha1, shape2 = beta1)*dbeta(y, shape1 = alpha2, shape2 = beta2), y,1)$value}) }, 0, 1)
This works for smaller alpha and beta even up to these values. But if I go much larger in beta the integrate function starts breaking. I've been comparing this to the Monte Carlo integration
n_sim = 1000000 # number of simulations y = rbeta(n_sim, shape1 = alpha2, shape2 = beta2) C2 = (1pbeta(y, shape1 = alpha1, shape2 = beta1)) mean(C2)
which yields approximately the same answers for smaller reasonable alpha and beta. And yes I need larger beta, I understand questioning why I need such large beta. It follows from large data sets. Does anybody know why this breaks? Maybe the inner workings of the integrate function? Is there a work around?

Integrate oscillatory functions with singularities
I am trying to evaluate the integral of the function
$\frac{1}{\sqrt{T}\sin(T/2)} \cos\left( \dfrac{\pi}{4} + \zeta*\left( k*T (y+\frac{T}{2})(1T/2\cot(T/2)) + \frac{T}{2}z+ \frac{x^2}{2T} + \frac{1}{4} \cot(T/2)(y^{2}+z^{2})\right) \right)$
in $(0,\infty)$ using (so far) scipy.integrate.quad. However, I always received the message:
IntegrationWarning: The integral is probably divergent, or slowly convergent.
At the moment, I have written in python:
def intreal(t): return 1/np.sqrt(t)/np.sin(t/2)*np.cos(np.pi/4 + zeta*(k*t  t/2*(1t/2*pow(np.tan(t/2),1)) + 1/(2*t)*x1*x1 + 1/4*pow(np.tan(t/2),1)*(pow(y1,2)+pow(z1,2))  (1t/2*pow(np.tan(t/2),1))*y1 + t/2*z1)) pts = [] for i in range (0,17): pts.append(2*i*np.pi) quad(intreal,1e12,33*np.pi,limit=500,points=pts,epsabs=0)
where intreal is the function defined above. I am also feeding the function with the values for x,y,z,k,\zeta.
In order to consider the singularities I have defined the problematic points, but I haven't had any luck at getting a reasonable result. Notice I would like to integrate up to infinity.

Quicker calculation of double integral in python (like MatLab's integral2)
I'm needing to perform a 2Dintegration (one dimension has an infinite bound). In MatLab, I have done it with integral2:
int_x = integral2(fun, 0, inf, 0, a, 'abstol', 0, 'reltol', 1e6);
In Python, I've tried scipy's dblquad:
int_x = scipy.integrate.dblquad(fun, 0, numpy.inf, lambda x: 0, lambda x: a, epsabs=0, epsrel=1e6)
and have also tried using nested single quads. Unfortunately, both of the scipy options take ~80x longer than MatLab's.
My question is: is there a different implementation of 2D integrals within Python that might be faster (I've tried "quadpy" without much benefit)? Alternatively, could I compile MatLab's integral2 function and call it from python without needing the MatLab runtime (and is that even kosher)?
Thanks in advance! Brad
Update:
Turns out that I don't have the "reputation" to post an image of the equation, so please bear with the formatting: fun(N,t) = P(N) N^2 S(N,t), where P(N) is a lognormal probability distribution and S(N,t) is fairly convoluted but is an exponential in its simplest form and a hypergeometric function (truncated series) in its most complex form. N is integrated from 0 to infinity and t from 0 to pi.