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

Howto compress a ellipsoid point cloud back to a sphere
I'm playing around with a ecompass chip (lsm303c). The chip returns vectors that, in a perfect setup, should represent a sphere around the origin. However due to hardiron and softiron effects the sphere is distorted. Hardiron distortion results in the sphere being displaced from the center, softiron distortion results in the sphere becoming a rotated ellipsoid. To compensate for these effects the chip has to be calibrated.
I found Matlab port of bounding ellipsoid and Bounding ellipsoid that describe a way to fit a pointset into a ellipsoid. It is based on Matlab code by code Nima Moshtagh here. A copy of the code from Matlab port of bounding ellipsoid is listed below.
Function
mvee
will returnA, c
.c
is the displacement vector andA
is the equation of the ellipse in matrix form.The question: I'm searching for the matrix that would transform the ellipsoid/ellipse shaped pointset that was used as input to
mvee
back to a sphere/circle like shape. How can I obtain that transformation matrix fromA
. It should not rotate but only compress the ellipsoid/ellipse pointset back to spheric/circle form. I guess this is nontrivial.import numpy as np import numpy.linalg as la import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D pi = np.pi sin = np.sin cos = np.cos def mvee(points, tol = 0.001): """ Finds the ellipse equation in "center form" (xc).T * A * (xc) = 1 """ N, d = points.shape Q = np.column_stack((points, np.ones(N))).T err = tol+1.0 u = np.ones(N)/N while err > tol: # assert u.sum() == 1 # invariant X = np.dot(np.dot(Q, np.diag(u)), Q.T) M = np.diag(np.dot(np.dot(Q.T, la.inv(X)), Q)) jdx = np.argmax(M) step_size = (M[jdx]d1.0)/((d+1)*(M[jdx]1.0)) new_u = (1step_size)*u new_u[jdx] += step_size err = la.norm(new_uu) u = new_u c = np.dot(u,points) A = la.inv(np.dot(np.dot(points.T, np.diag(u)), points)  np.multiply.outer(c,c))/d return A, c #some random points points = np.array([[ 0.53135758, 0.25818091, 0.32382715], [ 0.58368177, 0.3286576, 0.23854156,], [ 0.18741533, 0.03066228, 0.94294771], [ 0.65685862, 0.09220681, 0.60347573], [ 0.63137604, 0.22978685, 0.27479238], [ 0.59683195, 0.15111101, 0.40536606], [ 0.68646128, 0.0046802, 0.68407367], [ 0.62311759, 0.0101013, 0.75863324]]) # Singular matrix error! # points = np.eye(3) A, centroid = mvee(points) U, D, V = la.svd(A) rx, ry, rz = 1./np.sqrt(D) u, v = np.mgrid[0:2*pi:20j, pi/2:pi/2:10j] def ellipse(u,v): x = rx*cos(u)*cos(v) y = ry*sin(u)*cos(v) z = rz*sin(v) return x,y,z E = np.dstack(ellipse(u,v)) E = np.dot(E,V) + centroid x, y, z = np.rollaxis(E, axis = 1) fig = plt.figure() ax = fig.add_subplot(111, projection='3d') ax.plot_surface(x, y, z, cstride = 1, rstride = 1, alpha = 0.05) ax.scatter(points[:,0],points[:,1],points[:,2]) plt.show()

Why does ordering matter for Matlab?
I'm trying to define two variables, dr and dc, which store the dimensions of my matrix "data." It has 7 rows and 4 columns, so I expect to assign dr to 7 and dc to 4.
Why does
size(data) = [dr,dc]
Not work, while
[dr,dc] = size(data)
Work just fine?

How to go to the other round of two loops at the same time?
How to go to the other round of two loops at the same time?
Code:
parameters = {'CT', 'Imp', 'F1', 'F2'}; Time = {'T0', 'T1', 'T2', 'T3', 'T4', 'T5', 'T6'}; for i_parameters = 1: numel(parameters) my_parameters = parameters{i_parameters}; for i_Time = 1 : numel (Time) my_time = Time{i_Time}; for j = 1 : 7 Difference.(my_parameters).(my_time) = [Diff.(my_parameters)(:,j); Diff.(my_parameters)(:,j+7); Diff.(my_parameters)(:,j+14); Diff.(my_parameters)(:,j+21); Diff.(my_parameters)(:,j+28)]; DiffMean.(my_parameters).(my_time) = mean(Difference.(my_parameters).(my_time)); Diffstd.(my_parameters).(my_time) = std(Difference.(my_parameters).(my_time)); end end end
I would like to leave the loops
i
andi_Time
at each turn and at the same time that I would like to apply:
i_Time = T0
thenj = 1 / i_Time = T1
thenj = 2 / i_Time = T2
thenj = 3 / i_Time = T3
thenj = 4 / i_Time = T4
thenj = 5 / i_Time = T5
thenj = 6 / i_Time = T6
thenj = 7
. 
Finding approximate value of 'e' in C language
I am trying to find the approximate value of 'e' using C program. I am getting the value of e to be true no matter what I enter the value of n. Please let me find the error in my code.
#include <stdio.h> int f; int k; int factorial (int f); int main () { int n,i; int e = 1; printf("Enter the value of n:"); scanf("%d",&n); for (i = 1; i <= n; i++) { e = e + (1/factorial(i)); } printf("The value of e is %d",e); return(0); } int factorial (int f) { for (k = 1; k <= f; k++) { f = f*k; } return(p); }

julia: efficient and tidy way to avoid many function arguments
I have been writing stochastic PDE simulations in Julia, and as my problems have become more complicated, the number of independent parameters has increased. So what starts with,
myfun(N,M,dt,dx,a,b)
eventually becomes
myfun(N,M,dt,dx,a,b,c,d,e,f,g,h)
and it results in (1) messy code, (2) increased chance of error due to misplaced function arguments, (3) inability to generalise for use in other functions.
(3) is important, because I have made simple parallelisation of my code to evaluate many different runs of the PDEs. So I would like to convert my functions into a form:
myfun(args)
where args contains all the relevant arguments. The problem I am finding with Julia, is that creating a
struct
containing all my relevant parameters as attributes slows things down considerably. I think this is due to the continual accessing of the struct attributes. As a simple (ODE) working example,function example_fun(N,dt,a,b) V = zeros(N+1) U = 0 z = randn(N+1) for i=2:N+1 V[i] = V[i1]*(1dt)+U*dt U = U*(1dt/a)+b*sqrt(2*dt/a)*z[i] end return V end
If I try to rewrite this as,
function example_fun2(args) V = zeros(args.N+1) U = 0 z = randn(args.N+1) for i=2:args.N+1 V[i] = V[i1]*(1args.dt)+U*args.dt U = U*(1args.dt/args.a)+args.b*sqrt(2*args.dt/args.a)*z[i] end return V end
Then while the function call looks elegant, it is cumbersome to rework accessing every attribute from the class and also this continual accessing of attributes slows the simulation down. What is a better solution? Is there a way to simply 'unpack' the attributes of a struct so they do not have to be continually accessed? And if so, how would this be generalised?
edit: I am defining the struct I use as follows:
struct Args N::Int64 dt::Float64 a::Float64 b::Float64 end

Solving Multidimentional ODE in Python
I have to solve numerically ODE of 2nd order. The issue is that I need to solve it simultaneously for two objects (in general N objects) and each object has 4 variables (2 dimensions): 2 in position and 2 in velocity. Please see code below to see the example.
When trying to solve it in Python (Scipy) with provided methods I always encounter the error  "raise ValueError("
y0
must be 1dimensional.") ValueError:y0
must be 1dimensional."What is required for y0 is something like:
array([[[1, 0], [0, 1]], [[2, 0], [0, 2]]])
I think it might be possible to do so in Matlab, but I would like to stay with Python. Also, I want to solve it by RungeKutta model.
The possible solution is to reshape everything into 1 dimension, reshape it back and return reshaped again, but it is probably not that effective and efficient.
Thank you for your help. Please ask me if something is not clear

Integrating over a constant function
I am trying to integrate over a constant function in MATLAB 2017a, but I am stuck. First of all when I integrate using the following script, I get the right output. So the script works for a
x0
which depends ont
.function E=sol(n,k) x0 = @(t) t^(2); j = 0; E = zeros(n,1); while j < n+1 ; K = matlabFunction(subs(po(j,k))) ; eval(sprintf('x%d = integral(K,0,1);',j+1)) ; E(j+1,1) = subs(sprintf('x%d',j+1)) j = j+1; end end
Where function
po(j,k)
is as follows,function A_j = po(j,k) % Adomian polynomials if j >0 x = sym('x',[1 j]); syms p; % Assinging a symbolic variable for p syms x0; S = x0+ sum(p.^(1:j) .* x) ; % Sum of p*x up to order j Q =f(S,k); % Taking the kth power of S, i.e. A_nc = diff(Q,p,j)/factorial(j); % Taking the jth order derivative A_j = subs(A_nc,p,0) ; % Filling in p=0 else syms x0; S = x0; A_j =f(S,k); % Taking the kth power of S, end end
And where
f(x,k)
is,function F = f(x,k) % Nonlinear function of k power F = x^k ; end
Now when I cal
sol(n,k)
it does work. But when I try to change myx0
function insol(n,k)
in a constant function like,function E=solcon(n,k) x0 = @(t) 2.*ones(size(t)); j = 0; E = zeros(n,1); while j < n+1 ; K = matlabFunction(subs(po(j,k))) ; eval(sprintf('x%d = integral(K,0,1);',j+1)) ; E(j+1,1) = subs(sprintf('x%d',j+1)) j = j+1; end end
It does not work, As you can see I added
*ones(size(t));
just to make it a function oft
. But unfortunately it still doesn't work when I call,K = matlabFunction(subs(po(j,k))) ;
I get,
@()4.0
And then I get an error when I call,
eval(sprintf('x%d = integral(K,0,1);',j+1))
Could anyone help me out trying to integrate over a constant?
The error I get when I call
solcon(10,2)
isError using symengine>@()4.0 Too many input arguments. Error in integralCalc/iterateScalarValued (line 314) fx = FUN(t); Error in integralCalc/vadapt (line 132) [q,errbnd] = iterateScalarValued(u,tinterval,pathlen); Error in integralCalc (line 75) [q,errbnd] = vadapt(@AtoBInvTransform,interval); Error in integral (line 88) Q = integralCalc(fun,a,b,opstruct); Error in solcon1 (line 7) eval(sprintf('x%d = integral(K,0,1);',j+1)) ;
EDIT 2 I used the following script,
function E=solcon(n,k) x0 = @(t) 2.*ones(size(t)); j = 0; E = zeros(n,1); while j < n+1 ; K = matlabFunction(subs(po(j,k))) ; fstr= func2str(K) if fstr(3) == ')'; x{j+1} = K*(10) else x{j+1} = integral(K,0,1) end E(j+1,1) = subs(x{j+1},1); j = j+1 end end
But the following error occurs,
Undefined operator '*' for input arguments of type 'function_handle'. Error in solcone1 (line 9) x{j+1} = K*(10);

Integrating over a constant in Matlab
I am trying to integrate over a constant function in Matlab 2017a, but I am stuck. First of all when I integrate using the following script, I get the right output.
function E=sol(n,k) x0 = @(t) t^(2); j = 0; E = zeros(n,1); while j < n+1 ; K = matlabFunction(subs(po(j,k))) ; eval(sprintf('x%d = integral(K,0,1);',j+1)) ; E(j+1,1) = subs(sprintf('x%d',j+1)) j = j+1; end end
Where function po(j,k) is as follows,
function A_j = po(j,k) % Adomian polynomials if j >0 x = sym('x',[1 j]); syms p; % Assinging a symbolic variable for p syms x0; S = x0+ sum(p.^(1:j) .* x) ; % Sum of p*x up to order j Q =f(S,k); % Taking the kth power of S, i.e. A_nc = diff(Q,p,j)/factorial(j); % Taking the jth order derivative A_j = subs(A_nc,p,0) ; % Filling in p=0 else syms x0; S = x0; A_j =f(S,k); % Taking the kth power of S, end end
And where f(x,k) is,
function F = f(x,k) % Nonlinear function of k power F = x^k ; end
Now when I cal sol(n,k) it does work. But when I try to change my x0 function in sol(n,k) in a constant function like,
function E=solcon(n,k) x0 = @(t) 2.*ones(size(t)); j = 0; E = zeros(n,1); while j < n+1 ; K = matlabFunction(subs(po(j,k))) ; eval(sprintf('x%d = integral(K,0,1);',j+1)) ; E(j+1,1) = subs(sprintf('x%d',j+1)) j = j+1; end end
It does not work, As you can see I added *ones(size(t)); just to make it a function of t. But when I call,
K = matlabFunction(subs(po(j,k))) ;
I get,
@()4.0
And then I get an error when I call,
eval(sprintf('x%d = integral(K,0,1);',j+1))
Could anyone help me out trying to integrate over a constant ? with kind regards,
EDIT,
The error I get when I call solcon(10,2) is
Error using symengine>@()4.0 Too many input arguments. Error in integralCalc/iterateScalarValued (line 314) fx = FUN(t); Error in integralCalc/vadapt (line 132) [q,errbnd] = iterateScalarValued(u,tinterval,pathlen); Error in integralCalc (line 75) [q,errbnd] = vadapt(@AtoBInvTransform,interval); Error in integral (line 88) Q = integralCalc(fun,a,b,opstruct); Error in solcon1 (line 7) eval(sprintf('x%d = integral(K,0,1);',j+1)) ;

How does the gfortran compiler handle numerical coefficients in expressions?
I'm currently working on a library of temporal integrators. I'm wondering what's best for computing efficiency: Computing the coefficients for every time step, or using them as parameters from an external repository. The first case would be this way:
Subroutine getRKF45(this, dt, y, z) Implicit none Class(diffEc), Intent(InOut) :: this Real(8), Intent(In) :: dt Real(8), dimension(:), Intent(InOut) :: y, z Real(8), dimension(6,this%n) :: k Call field(this%deriv, this%x, this%t) k(1,:) = this%deriv(:) Call field(this%deriv, this%x+k(1,:)*dt*(1.d0/5.d0), this%t+dt*(1.d0/5.d0)) k(2,:) = this%deriv(:) Call field(this%deriv, this%x+k(1,:)*dt*(3.d0/40.d0)+k(2,:)*dt*(9.d0/40.d0) & , this%t + dt*(3.d0/10.d0)) k(3,:) = this%deriv(:) Call field(this%deriv, this%x+k(1,:)*dt*(3.d0/10.d0)k(2,:)*dt*(9.d0/10.d0) & +k(3,:)*dt*(6.d0/5.d0), this%t + dt*(3.d0/5.d0)) k(4,:) = this%deriv(:) Call field(this%deriv, this%xk(1,:)*dt*(11.d0/54.d0)+k(2,:)*dt*(5.d0/2.d0) & k(3,:)*dt*(70.d0/27.d0)+k(4,:)*dt*(35.d0/27.d0), this%t + dt) k(5,:) = this%deriv(:) Call field(this%deriv, this%x+k(1,:)*dt*(1631.d0/55296.d0)+k(2,:)*dt*(175.d0/512.d0) & +k(3,:)*dt*(575.d0/13824.d0)+k(4,:)*dt*(44275.d0/110592.d0) & +k(5,:)*dt*(253.d0/4096.d0), this%t + dt*(7.d0/8.d0)) k(6,:) = this%deriv(:) y(:) = this%x(:) + dt*(k(1,:)*(37.d0/378.d0)+k(3,:)*(250.d0/621.d0)+k(4,:)*(125.d0/594.d0) & +k(6,:)*(512.d0/1771.d0)) z(:) = this%x(:) + dt*(k(1,:)*(2825.d0/27648.d0)+k(3,:)*(18575.d0/48384.d0) & +k(4,:)*(13525.d0/55296.d0)+k(5,:)*(277.d0/14336.d0)+k(6,:)*(1.d0/4.d0)) End Subroutine getRKF45
And using coefficients precomputed:
Subroutine getRKF45(this, dt, y, z) Use IntegratorUtilities, only: rk45Coef Implicit none Class(diffEc), Intent(InOut) :: this Real(8), Intent(In) :: dt Real(8), dimension(:), Intent(InOut) :: y, z Real(8), dimension(6,this%n) :: k Call field(this%deriv, this%x, this%t) k(1,:) = this%deriv(:) Call field(this%deriv, this%x+k(1,:)*dt*rk45Coef(1), this%t+dt*rk45Coef(1)) k(2,:) = this%deriv(:) Call field(this%deriv, this%x+k(1,:)*dt*rk45Coef(2)+k(2,:)*dt*rk45Coef(3) & , this%t + dt*rk45Coef(4)) k(3,:) = this%deriv(:) Call field(this%deriv, this%x+k(1,:)*dt*rk45Coef(4)k(2,:)*dt*rk45Coef(5) & +k(3,:)*dt*rk45Coef(6), this%t + dt*rk45Coef(7)) k(4,:) = this%deriv(:) Call field(this%deriv, this%xk(1,:)*dt*rk45Coef(8)+k(2,:)*dt*rk45Coef(9) & k(3,:)*dt*rk45Coef(10)+k(4,:)*dt*rk45Coef(11), this%t + dt) k(5,:) = this%deriv(:) Call field(this%deriv, this%x+k(1,:)*dt*rk45Coef(12)+k(2,:)*dt*rk45Coef(13) & +k(3,:)*dt*rk45Coef(14)+k(4,:)*dt*rk45Coef(15) & +k(5,:)*dt*rk45Coef(16), this%t + dt*rk45Coef(17)) k(6,:) = this%deriv(:) y(:) = this%x(:) + dt*(k(1,:)*rk45Coef(18)+k(3,:)*rk45Coef(19)+k(4,:)*rk45Coef(20) & +k(6,:)*rk45Coef(21)) z(:) = this%x(:) + dt*(k(1,:)*rk45Coef(22)+k(3,:)*rk45Coef(23) & +k(4,:)*rk45Coef(24)+k(5,:)*rk45Coef(25)+k(6,:)*rk45Coef(26)) End Subroutine getRKF45
I've run some test and the time difference between the two schemes is not relevant. Does the compiler preprocess these coefficients multiplications?