# Numerical Integration [duplicate]

Let $\alpha=1.714$ and $\lambda=-v(1-v^{\alpha})/\alpha$, where $0<v<1$. Consider $\phi(v)=\exp\{\int_{v_{0}}^{v}\frac{1}{\lambda(t)} dt \}$, where $0< v_{0} <1$ is an arbitrary chosen constant. How to integrate the given integral numerically and plot $\phi(v)$ vs $v$ using Matlab?

• Welcome to Scicomp.SE! I think this has already been answered here: scicomp.stackexchange.com/questions/21366/…. You should look at that question, and if the answers do not help you, edit your question to describe what precisely you need beyond that. – Christian Clason Aug 23 '16 at 10:19

The integral function in MATLAB performs numerical integration. See the MATLAB help, it's straightforward to use.

Actually, an analytical solution of your integral exists. The calculations are a bit tricky but feasible. Let's rewrite your integral in your exponential (assuming $v=v(t)$) :

$$I = \alpha \int_{v_0}^v{\frac{dt}{v(v^{\alpha}-1)}}$$

We perform a change of variable : $$u=v^{\alpha}-1$$ $$du = \alpha v^{\alpha-1}dt$$ $$u_0 = v_0^{\alpha}-1$$ then $I$ becomes : $$I = \alpha \int_{u_0}^u{\frac{du}{\alpha~v^{\alpha-1}~v~u}} = \int_{u_0}^u\frac{du}{v^{\alpha}~u} = \int_{u_0}^u\frac{du}{(u+1)u}$$ Much better to integrate right ?

Now an element decomposition gives you : $$I = \int_{u_0}^u\frac{du}{u} - \int_{u_0}^u\frac{du}{u+1}$$ And the final blow: $$I = [ln(u)]_{u_0}^u - [ln(u+1)]_{u_0}^u$$ Now going back to $v$ : $$I = ln \left(\frac{v^\alpha-1}{v_0^\alpha-1}\left(\frac{v_0}{v}\right)^\alpha \right)$$ So taking the exponential gives you the analytical solution of $\phi(v)$:

$$\phi(v) = \frac{v^\alpha-1}{v_0^\alpha-1}\left(\frac{v_0}{v}\right)^\alpha$$

The code below shows the comparison between your numerical integration and the analytical solution.

clc;
clear all
close all;

v=[0.1:0.01:0.9]; % Vector v
v0=0.05; % Initial condition
a=1.714; % alpha

fun = @(v) 1./(-v.*(1-v.^a)./a); % lambda

for i=1:length(v)
phi_num(i) = integral(fun,v0,v(i)); % Numerical integration
end

phi_ex=((v0./v).^a) .* ((v.^a-1)./(v0.^a-1)); % Analytical solution

plot(v,exp(phi_num),'b.')
hold on
plot(v,phi_ex,'r')
legend('Numerical','Analytical')