How to integrate this function numerically?

Hello everyone,

I'm stuck at this equation and I want to integrate it numerically

where

and r = 695510, a = -0.0021, b = 1.34, Wo = 438.1

I need to integrate the first equation numerically to get V as a function of R

Numerical integration from R = 10, where it is assumed V = U, to 215 should give V(R).

I appreciate your help.

Thank you,

[Walter Roberson]

Is e the base of the natural logarithms? And is it acting as a constant multiplier or is e being raised to the part after it?

[Mohamed Nedal]

Yes that's right, it's the base of the natural logarithm and it's being raised to the part after it.

Answer 1

[Chuguang Pan]

You can use Euler formula. V(n+1)=V(n)+h*f(R,V(n)), which f(R,V) is the right side of differential equation.

But you need to know the initial value V(10).

r=695510;a=-0.0021;b=1.34;Wo=438.1;
h=0.1;%Integral step size, you can change this value
N=100;%It means that you want to integral from 10 to 10+h*N
V=zeros(1,N+1);%Initialization V array
V(1)=?;%Need to know the initial value V(10)
for n=1:N
    R=10+(n-1)*h;
    W=Wo*sqrt(1-exp(1)*(2.8-R)/8.1);
    V(n+1)=V(n)+h*r*a*R^(-b)*(1-W/V(n));
end
X=10+(0:N)*h;
Y=V;
plot(X,Y);

[Mohamed Nedal]

Thanks for your response. The initial value of V(10) is known (693).

I think the 3rd line of the for loop should be like this

W=Wo*sqrt(1-exp((2.8-R)/8.1)); 

Could you please elaporate on h and N?

What are the consequences of changing the integral step size (h) and N?

[Chuguang Pan]

r=695510;a=-0.0021;b=1.34;Wo=438.1;
h=0.1;%Integral step size, you can change this value
N=(215-10)/h;%It means that you want to integral from 10 to 215
V=zeros(1,N+1);%Initialization V array
V(1)=693;%Need to know the initial value V(10)
for n=1:N
    R=10+(n-1)*h;
    W=Wo*sqrt(1-exp((2.8-R)/8.1));
    V(n+1)=V(n)+h*r*a*R^(-b)*(1-W/V(n));
end
X=10+(0:N)*h;
Y=V;
plot(X,Y);

you can decrease h for increasing presision

[Mohamed Nedal]

Thank you so much, that's very helpful.

[Mohamed Nedal]

Please I have another question. Knowing the above parameters, How do we get t from the folowing equation?

[Chuguang Pan]

If t is independent variable, it is similar to R. You need to decide for youself.

Such as t=t0+(0:N)*h.

[Mohamed Nedal]

t is the time (in seconds) the object takes to cross the distance R at speed V.

to = 0, I need to find t(R=215).

[Mohamed Nedal]

Could you check this please.

r = 695510; a = -0.0021; b = 1.34; Wo = 438.1; 
h = 1;              % Integral step size, you can change this value 
N = (215 - 10)/h;   % It means that you want to integral from 10 to 215 
V = zeros(1,N+1);   % Initialization V array 
t = zeros(1,N+1); 
t(1) = 0; 
V(1) = 693;         % Need to know the initial value V(10) 
for n = 1:N 
    R = 10 + (n-1)*h; 
    W = Wo * sqrt(1-exp((2.8-R)/8.1)); 
    V(n+1) = V(n) + h*r*a*R^(-b)*(1-W/V(n)); 
    t(n+1) = t(n) + h*a*R^(-b)*(V(n)-W); 
end 
X = 10 + (0:N)*h; 
t = (0:N)*h;        % to = 0 
plot(X,V); 

t here is monotonically increasing (t=0,1,2,...) and when I removed this line

t = (0:N)*h; 

it gives negative fractions and it's not correct.

[Chuguang Pan]

V=dR/dt. dt=dR/V.

R=10+(0:N)*h;
dR=diff(R);
V=(V(2:end)+V(1:end-1))*0.5;%average speed
dt=dR./V;
t=[0 cumsum(dt)];

[Mohamed Nedal]

Thank you very much :)

[Mohamed Nedal]

Excuse me another question..

I would like to apply this equation on several objects and I tried to replace V(1), which's the initial value V(10), with a vector column instead of a single value like this one for instance

Vtest = [136
490
359
361
693
1374
585
542
1118
910
1188]; 

so that each value is assigned for one object. But I got this error

In an assignment  A(:) = B, the number of elements
in A and B must be the same. 

Error in numeric_integ (line 22)
V(1) = Vtest;         % Need to know the initial value V(10) 

Could you help me out?

and if I want to find only the value of V(215) and t(215), which are the last values of V and t at distance 215, how can I do that?

[Chuguang Pan]

result_V=zeros(1,length(Vtest));% save last V for every Vtest
result_t=zeros(1,length(Vtest));% save last t for every Vtest 
for i=1:length(Vtest) %for every Vtest
    V=zeros(1,N+1);
    V(1)=Vtest(i);
    %do the computing
    result_V=V(end); %using end to index last value of array V
    result_t=t(end); %using end to index last value of array t
end

[Mohamed Nedal]

Thank you so much