Path: news.mathworks.com!newsfeed-00.mathworks.com!news.tele.dk!feed118.news.tele.dk!news.tele.dk!small.news.tele.dk!newsfeed.xs4all.nl!newsfeed6.news.xs4all.nl!xs4all!feeder.news-service.com!94.75.214.39.MISMATCH!aioe.org!.POSTED!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: ode with definite integral
Date: Thu, 30 Jun 2011 01:41:00 -0700
Organization: Aioe.org NNTP Server
Lines: 104
Message-ID: <iuhcr1$jqq$1@speranza.aioe.org>
References: <iuh6ma$ggi$1@newscl01ah.mathworks.com> <iuh9ga$mua$1@newscl01ah.mathworks.com>
Reply-To: <HIDDEN>
NNTP-Posting-Host: TUXTYYqX1yG7hs3zxUg7ng.user.speranza.aioe.org
Mime-Version: 1.0
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 7bit
X-Complaints-To: abuse@aioe.org
User-Agent: Mozilla/5.0 (Windows; U; Windows NT 6.1; en-US; rv:1.9.2.18) Gecko/20110616 Thunderbird/3.1.11
X-Notice: Filtered by postfilter v. 0.8.2
Xref: news.mathworks.com comp.soft-sys.matlab:734293

On 6/30/2011 12:44 AM, Roger Stafford wrote:
> "Grzegorz Knor" wrote in message<iuh6ma$ggi$1@newscl01ah.mathworks.com>...
>> Hello,
>> I wonder how to solve such a task:
>> Suppose we want to solve a differential equation in the form:
>> dy/dt = f(t) = exp(-te)
>> where:
>> te = \int from 0 to t (2^(y/10)) dt
>> Is it possible to do it with Matlab ode solvers?
>>
>> Best regards
>> Grzegorz
> - - - - - - - - -
>    Yes, that can easily be done with an ode solver.  Define
>
>   te(t) = int('2^(y(s)/10)','s',0,t)
>
> Then we have
>
>   dte(t)/dt = 2^(y(t)/10)
>
> and
>
>   dy(t)/dt = exp(-te(t))
>
> so this is in the proper form to be solved with an ode solver.
>  You must of course arrange that te(0) = 0.  The initial value, y(0), is up to you.
>
> Roger Stafford

Nice Roger, that is a nice way to do it. Below I show your method,
implemented in ode45, using initial value for te(0)=0, and then
show another method I just hacked using the original formula as
is by using quad directly, and interp1 to do the integration
at each step. It is a hack, but it gives the same answer.

Had to protect against 2 special cases: when ode45 uses
the same time step twice, else intepr1 will not be happy,
and for initial case at t=0.

Here we go:

Roger method: (the right way to do it)

------------------
functionfoo1

figure;
y0=1;
integralAtZero=0;
maxTime = 10;
    
[t,y]=ode45(@RHS,[0 maxTime],[y0  integralAtZero]);
plot(t,y(:,1));
title('solution y(t)');

     function ydot = RHS(t,y)
         ydot = [exp(-y(2)) ; 2^(y(1)/10) ];
     end
     
end
---------------------------

my method, just for fun :)

------------------------------
function foo

y0=1;
maxTime = 10;
globalY = zeros(1000,1);    % change as needed
globalt = zeros(1000,1);   % to prevent dynamic allocation
globalydot = zeros(1000,1); %just  a hack
indx=0;
    
[t,y]=ode45(@RHS,[0 maxTime],y0);
plot(t,y);


     function ydot = RHS(currentTime,y)
         if currentTime(1)>0 && currentTime(1) == globalt(indx)
                ydot= globalydot(indx);
         else
           indx=indx+1;
           globalY(indx)=y(1);
           globalt(indx)=currentTime(1);
           ydot = exp(-quad(@f,0,currentTime));
           globalydot(indx)=ydot;
         end
     end

     function v = f(t) % integrate 2^(y(t)/10)
         if all(t==0)
            v=t;
         else
           yi = interp1(globalt(1:indx), globalY(1:indx),t) ;
           v = 2.^(yi/10);
         end
     end
end
-------------------------------------


--Nasser