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

change the plot in matlab
When a figure is created in MATLAB , Is it possible to change that shape?
Can I make such changes in the figure window? For example, I want to change the position of some nodes, or change the color of the arrows. Some angles are large, for example, in my shape
the angle between 5, 7 and 12 is enlarged and I want to be smaller.

Why does Fast tif loading cause a Dimension mismatch? MATLAB
I'm trying to debug this set of code but I do not understand why an error is occurring. I found a set of code released back in 2012 which was designed for fast loading of tiff stacks into matlab. However for some reason when I run the code I always receive a dimension mismatch error.
code in question. To get to run, need the
tifflib.mexw64
file in the directory path.FileTif='ImageStack.tif'; InfoImage=imfinfo(FileTif); mImage=InfoImage(1).Width; nImage=InfoImage(1).Height; NumberImages=length(InfoImage); FinalImage=zeros(nImage,mImage,NumberImages,'uint16'); FileID = tifflib('open',FileTif,'r'); rps = tifflib('getField',FileID,Tiff.TagID.RowsPerStrip); for i=1:NumberImages tifflib('setDirectory',FileID,i); % Go through each strip of data. rps = min(rps,nImage); for r = 1:rps:nImage row_inds = r:min(nImage,r+rps1); stripNum = tifflib('computeStrip',FileID,r); FinalImage(row_inds,:,i) = tifflib('readEncodedStrip',FileID,stripNum); end end tifflib('close',FileID);
I tried to look through for loops to see how the different loops are different but they all seem to be similar. And the code usually fails at the end of the first frame but the later cycle isn't any different from the first few.

Proportional Control Matlab
I am working on a boiler project in matlab with a partner. We are still new to this program and are currently working on the proportional control part of the boiler program. As a result we get 3 graphs like we are supposed to but we keep getting a mysterious dip on our center graph as seen below. Any diagnosis to the problem on how we can get rid of the dip is what we have been trying to figure out.
Here is our code:
clear all; close all; endtime=60; time=0; ResPer=0.05; %declares basic functions setpoint=7; flowrate=2.4; valveposition=0.2; delay=5; timeaxis=[1:endtime]; valveaxis=0; %declares arrays related to time and graphing steamaxis=0; propconstant=0.25; RelError=0; %declares arrays regarding the proportional controller desiredvalve=0.5823; while time<endtime %RelError=SteamRelativeError(flowrate, setpoint); %calculates the Relative Error used in the rest of the function if (time <= delay) valveposition=0.75*valveposition+0.25*(desiredvalvepropconstant*RelError); %deals with the valve position before the time exceeds the delay %valveaxis(time)=valveposition; else valveaxis(timedelay)=0.75*valveposition+0.25*(desiredvalvepropconstant*RelError); %calculates and stores valve position after the delay has been passed end RelError=SteamRelativeError(flowrate, setpoint); %calculates the Relative Error used in the rest of the function time=time+1; %increments time flowrate=ValvePerToFlowRate(valveposition); %calculates flow rate from valve position valveaxis(time)=valveposition; %stores temporary valve position steamaxis(time)=flowrate; %stores steam flow rate erroraxis(time)=RelError; %stores the Relative Error end subplot(3,1,2), plot(timeaxis, valveaxis, 'b'); axis([0 endtime 0.2 0.8]); xlabel('Time(min)'); ylabel('Valve%'); subplot(3,1,1), plot(timeaxis, steamaxis, 'b'); axis([0 endtime 2 setpoint+2]); xlabel('Time(min)'); ylabel('Steam mass flow rate (lbm/s)'); title('Proportional Constant=0.25'); %graphs functions subplot(3,1,3), plot(timeaxis, erroraxis, 'b'); axis([0 endtime 0.6 0.2]); xlabel('Time(min)'); ylabel('Error');
And here are our graphs:

Numerical solution to differential equation with multiple variables
dm4(t)/dt = γm3(t) − ρm4(t) where gamma and rho are known. with t_0=0 how would I use a numerical method to find some solutions to this differential equation?

Interpretation of the input p  Function gradient(f,p) in MATLAB
I'd like to know what is the input p in the function gradient(f,p). In particular, if f is a vector with n elements, each of them is a sample, and if the sample rate is sr=0.1 Hz, thus is p=1/sr=sample time?
Thank you for your time.

Riding the wave Numerical schemes for hyperbolic PDEs, lorena barba lessons, assistance needed
I am a beginner python user who is trying to get a feel for computer science, I've been learning how to use it by studying concepts/subjects I'm already familiar with, such as Computation Fluid Mechanics & Finite Element Analysis. I got my degree in mechanical engineering, so not much CS background.
I'm studying a series by Lorena Barba on jupyter notebook viewer, Practical Numerical Methods, and i'm looking for some help, hopefully someone familiar with the subjects of CFD & FEA in general.
if you click on the link below and go to the following output line, you'll find what i have below. Really confused on this block of code operated within the function that is defined.
Anyway. If there is anyone out there, with any suggestions on how to tackle learning python, HELP
In[9]
rho_hist = [rho0.copy()] rho = rho0.copy() **# im confused by the role of this variable here** for n in range(nt): # Compute the flux. F = flux(rho, *args) # Advance in time using LaxFriedrichs scheme. rho[1:1] = (0.5 * (rho[:2] + rho[2:])  dt / (2.0 * dx) * (F[2:]  F[:2])) # Set the value at the first location. rho[0] = bc_values[0] # Set the value at the last location. rho[1] = bc_values[1] # Record the timestep solution. rho_hist.append(rho.copy()) return rho_hist

Double integral in RcppNumerical
I have to calculate the double integral of the function:
> DIntegral < function(x,y){res < pnorm(x,1,0.1) * dexp(y2,1.2) return(res) }
The upper limit for x and y are:
10
andInfinity
respectively. The lower limit for x and y are:1
and2
respectively.How can I do this double integration in
RcppNumerical
?For one dimensional integration my C++ file looks like:
// [[Rcpp::depends(RcppEigen)]] // [[Rcpp::depends(RcppNumerical)]] #include <RcppNumerical.h> using namespace Numer; class PDF: public Func { private: double beta; double M0; public: PDF( double beta_, double M0_): beta(beta_), M0(M0_) {}; double operator()(const double& x) const { return R::dexp(xM0,beta,0); } }; // [[Rcpp::export]] double integrate_test2( double beta, double M0, double upper, double lower) { PDF f( beta, M0); double err_est; int err_code; double res = integrate ( f, lower, upper,err_est,err_code); return(res); }
Code for 2 dimensional integration with finite limits
// [[Rcpp::depends(RcppEigen)]] // [[Rcpp::depends(RcppNumerical)]] #include <RcppNumerical.h> using namespace Numer; class PDF: public MFunc { private: double mu; double sigma; double beta; double M0; public: PDF( double mu_, double sigma_, double beta_, double M0_): mu(mu_), sigma(sigma_), beta(beta_), M0(M0_) {}; double operator()(Constvec& x) { return R::pnorm(x[0],mu,sigma,1,0) * R::dexp(x[1]M0,beta,0); } }; // [[Rcpp::export]] double integrate_test2( double mu, double sigma, double beta, double M0) { Eigen::VectorXd lower(2); lower << 1, 2; Eigen::VectorXd upper(2); upper << 10, 50; PDF f( mu, sigma, beta, M0); double err_est; int err_code; double err_est2; int err_code2; double res = integrate ( f, lower, upper,err_est,err_code); return(res); }

Using numeric integration to find the area enclosed by loop of coordinates
I am trying to find or write some code that can give me the area between two lines, and if the lower one crosses the higher one an negative area is accounted for. For this (cannot post an image :( ) . I have writen code but am not 100% confident with the answers. I have also tried to just using a numpy function of np.trapz.
def TrapzCn(n): Uy = Upper(n) Ly = Lower(n) xc = Open(L[0]) Top = np.trapz(Uy,x=xc[:10]) Bottom = np.trapz(Ly,x=xc[:10]) return(Top + Bottom)
and here's mine
def Cn(n): Uy = Upper(n) xc = Open(L[0]) Ly = Lower(n) for i in range(len(Ly)): Ly[i] = float(Ly[i]) Uy[i] = float(Uy[i]) A = 0 i = 0 while i+1 < len(Uy): A1 = (xc[i+1]  xc[i]) * ((Uy[i] + Uy[i+1])/2) A2 = (xc[i+1]  xc[i]) * ((Ly[i] + Ly[i+1])/2) if A1 > A2 and Uy[i] < 0.0 : A += abs(A1) + abs(A2) elif A1 > A2 and ((Uy[i] + Uy[i+1] )/2)> 0: A += abs(A1)  abs(A2) elif A1 < A2 and A1 < 0.0 : A = (abs(A2)  abs(A1)) i += 1 return(A)
Where the Upper and Lower functions are just unpacking my data. So my main question is, does the np.trapz function account for the curve doubling back on itself as well as the negative area once they cross. And if not, Does mine?

Write a trapezoidal rule routine to integrate the following Bessel integral using MATLAB without builtin functions
I am trying to create a trapezoidal rule routine to integrate the following Bessel integral in Matlab without using builtin functions. I have only solved integrals with one variable, but now there are two variables. I am not sure how to approach the problem provided in the link below. This is a screenshot of the problem