solve integral using ode45 (by derivative)

Hi!

I've got a task that requires me to solve the following integral using ode45 by derivating the integral with the upper limit as "t". I don't seem to understand how to solve this, other than using the command "integral" (which I'm not allowed to do).

I've been given m_e, but my ODE45 attempt doesn't seem to work, and I don't understand why. The code below is what I've come up with.

planckfunnew=@(t,T)planck(lambda,T);
[tsome,ssome]=ode45(planckfunnew,[0 1],[4.*10.^-7;100]);
% plot(ssome(:,1),ssome(:,2))

function me=planck(lambda,T)
h=6.6256e-34;
c=2.9979e8;
k=1.3805e-23;
a=2*pi*h*c^2;
b=h*c/k./lambda./T;
me=a./lambda.^5./(exp(b)-1);

I'm not that sure of how to define the function, which I'm supposed to insert in the ode45 line.

I'm grateful for any help that I recieve!

Answer 1

[John D'Errico]

While what Ameer wrote looks like it will work, it still uses the dreaded & disallowed integral function. I'm adding this answer because I think many people can benefit from the ideas of how to use ODE45 as a general integration tool, to explain what is needed. In fact, I can think of a few cases where I have indeed used ODE45 to do a numerical integration when others might have just looked at integral. As well, I would add that what you tried would definitely not work, but it does move in the general desired direction.

The trick is to see that ODE45 is just a variation of a numerical integration code, that integrates starting from one end of the interval (typically the left) and works to the right. tspan gives the limits of integration.

Here lambda is the independent variable, and T is a parameter that we need to make the result a function thereof. Be careful, as it is easy to get variable names confused, even though ODE45 uses t in the call. Here lambda takes the place of t in ODE45. The clue there is the d(lambda) in the integral.

You will be integrating M_e(lambda,T), with M_s(T) being the final result. And apparently, we are not asked to see the integral as a function of lambda, but only the final value at the end of the interval. With all of that explained, let me now explain the code I will write.

The function handle m_e encapsulates the current values of each of h, c, k, and a. However it is a function of lambda and T.

What does ODE45 expect though? It wants something of the general form

dy/dt = F(t,y)

Again, BE CAREFUL! Here t is lambda, but y is the unknown, M_s. (I think that was an s. The picture of your equation was just a bit fuzzy though. Hard on these bleery eyes.) In fact, M_s never appears in the kernel of the integration, so we need not worry about it.

One other problem. An ODE solver, just like any integration includes an arbitrary constant of integration. When you supply the initial value for the ODE solver, that allows it to resolve the arbitrary constant. It seems logical to call that value 0 here, but I'm not the person to know, and I'm feeling far too lazy to do some real thinking for the moment. (Hey, its hot out today!) With all of that explained...

function M_s = myplanck(T)
  % Some constants
  h=6.6256e-34;
  c=2.9979e8;
  k=1.3805e-23;
  a=2*pi*h*c^2;
  
  % easiest to just write m_e in one line as a function handle
  m_e = @(lambda,T) a./lambda.^5./(exp(h*c/k./lambda./T)-1);
  
  % we are almost ready to call ODE45. Set up the tspan (ok, lambdaspan)
  % and the initial condition. Then create the function as needed by ODE45.
  M_s0 = 0;
  tspan = [4e-7,7e-7];
  ODEFUN = @(lambda,y) m_e(lambda,T);
  
  % The tilde as the first output argument means we don't care what tout is,
  % because we are not interested in the integral as a function of lambda,
  % only the full integral over tspan from start to finish.
  [~,yout] = ode45(ODEFUN,tspan,M_s0);
  
  % M_s is the final element of yout, thus the integral over the entire interval.
  M_s = yout(end);
end

Now what are the odds I did this all, while making no overt mistakes? Sometimes I get lucky. Really though, luck favors the prepared. And when you write code very carefully, it can work the first time. You need to think it through to the end, and visualize how everything interacts. Then the entire thing becomes easy.

myplanck(10000)
ans =
   1.8551e+08

That seems consistent with what Ameer wrote, using integral.

The trick is to understand what the ODE solver wants, as well as to see through the variable names. Don't allow yourself to get confused by things. You also need to see how the variables get passed around, for example, what uses T, and how does that variable get passed into the necessary functions.

Finally, be careful. Suppose you next wanted to make this code work for multiple values of T? While the code I wrote will fail if you passed in a vector of values for T, the simplest solution would be to create a loop inside the myplanck function. The loop would operate on each element of the vector T, one at a time. So some slight additional complexity, but not a lot.

Answer 2

[Brent Kostich]

It appears you are working with the Stephan-Boltzmann Law. However, I think the syntax in your "planck" function may be a bit off. In any case, here is what I think you are trying to accomplish. The main function accepts a range of temperature values, used to evaluate your given equations.

function [] = wavelength_fractal(T)

%{
  T should typically be a vector input of temperature range
%}

%% Calculating Spectral Intensity
% pre-allocating a vector of length T
M_s = zeros(1, length(T));

for i = 1:length(T)
    % numerical integration "planck" function
    [~, y] = ode45(@(lambda, y) planck(lambda, T(i)), [4e-7 7e-7], 0);
    
    M_s(i) = y(end);
end

%% Using Stefan-Boltzmann constant
stef_bolt = 5.67e-8;
M_e = stef_bolt*T(i)^4;

%% Calculating fraction
wave_frac = M_s/M_e;

%% Plotting Output
figure(1)

plot(T, wave_frac)
grid

function m_e = planck(lambda, T)
h = 6.6256e-34;
c = 2.9979e8;
k = 1.3805e-23;

a = 2*pi*h*c^2;
b = (h*c)/(k*lambda*T);

m_e = a/(lambda^5*(exp(b)-1));

So then by entering a command line prompt of:

wavelength_fractal(linspace(100,10000))

The result should be:

You of course can add your own figure title, axis labels, etc. to spruce it up. Cheers!

[John D'Errico]

Thanks. I'll completely agree that I have no idea what the real equations needed here are, only that I worked from what the OP provided.

[Brent Kostich]

I hear you, your solution was completely sound to that end. At one point this thread contained more information about the OP's problem, but I think it has since been taken down. My solution was considering that additional information.

Answer 3

[Ameer Hamza]

Try this

M = @(T) integral(@(lambda) planck(lambda,T), 4e-7, 7e-7);
s = M(10000); % T=10000

function me=planck(lambda,T)
h=6.6256e-34;
c=2.9979e8;
k=1.3805e-23;
a=2*pi*h*c^2;
b=h*c/k./lambda./T;
me=a./lambda.^5./(exp(b)-1);
end

Result

>> s
ssome =
   1.8551e+08