MATLAB - Using solution of one ODE in another

[RedAnsar]

Hi all,

I'm trying to use MATLAB to obtain simulations for some equations that describe a model. I'm new to MATLAB (though I've taken a course in C++ and another in Java), so I read a bit on the mathworks website on solving ODEs, and settled on ode45

The equations I'm trying to model are the following (see link):
http://i.imgur.com/vWATx.jpg

As can be seen, the first equation is a straightforward (I guess?) ODE, where the L function is defined, for some x, as L(x)=coth(x)-(1/x). I was not sure how to write that (d ln L/dE)^(-1), so I just had WolframAlpha evaluate it and I typed it in.
Mentioned earlier, I believe that I got the code down for the first ODE:

%this is the M-file firstode
%tau_B = 0.0016
function Eprime = firstode(t,E);
Eprime=(((0.0016).^(-1))*((coth(E)-(1./E))/((E).^(-2)-(csch(E)).^(2)))*(((10)/E)*cos(625*t)-1));

%this is the main file, firstodesolver.m
tspan=[0,0.2];
E0=10;
[t,E]=ode45(@firstode,tspan,E0);
plot(t,E, 'b-')

Assuming this gives me the right value for E (I didn't know how to write xi_E on MATLAB), how do I go about to using the solution for E in the second equation? I thought defining the solution, i.e. something like
sol = ode45(@firstode,tspan,E0)
would do it...but that doesn't seem to work. Can y'all provide any suggestions?
NOTE: The code written does not include anything re. the second equation as I've yet to solve the dilemma I'm facing.

Sorry for the long post. Thanks all!
MATLAB's a bit strange -- this would be nicer on C++ :(RedAnsar

 

[chiro]

Hey RedAnsar and welcome to the forums.

One suggestion I have is if you can solve a system of ODE's, then create a new function g(L) = ln(L) and then replace the d(LnL)/de with dg/de.

 

[RedAnsar]

Thanks for the welcome, chiro.

I thought it would be a system of ODEs, but it seems that would be only in the case that my equations wouldn't contain any terms with t, only the variable sought after
E.g.)
y1' = y1 + y2
y2' = y1*y2
In my case, I do have terms with t (notice the cos(w*t) in the equations)

Re. declaring that function...so, in C++, one can simply declare a variable a certain type
E.g.)

double x;
int y;

but can you do that as well in MATLAB?? I wanted to define L(x)=coth(x)-(1./x) but MATLAB tells me I need to define x. However, I don't want to give x an explicit range of values, just sort of say "x can be any number you plug in." Is that possible?

Again, thanks and sorry for lots of (probably trivial) questions

Red Ansar

 

[chiro]

The main data-type in MATLAB is a matrix and the matrix has to have values in it and not variables as far as I understand, and as far as my experience has told me.

The kind of thing you are talking about is dealt with computational math packages that use symbolic engines and these include things like Maple and Mathematica.

In terms of solving this thing in MATLAB, what I would recommend is that you implement your own numeric solver routine based on Runge-Kutta and in MATLAB depending on your delta-apprimoximation (i.e. deltat to approximate dt), you will generate a vector of values for the t's with a step-size of delta t and create another vector for the function values that are evaluated with runge-kutta scheme.

 

[RedAnsar]

In terms of solving this thing in MATLAB, what I would recommend is that you implement your own numeric solver routine based on Runge-Kutta and in MATLAB depending on your delta-apprimoximation (i.e. deltat to approximate dt), you will generate a vector of values for the t's with a step-size of delta t and create another vector for the function values that are evaluated with runge-kutta scheme.

Do you have a suggestion on how to go about doing this...? Sounds...daunting : (
I suppose a first step is reading more about Runge-Kutta, huh (it was never covered in my ODE class)

 

[chiro]

Do you have a suggestion on how to go about doing this...? Sounds...daunting : (
I suppose a first step is reading more about Runge-Kutta, huh (it was never covered in my ODE class)

Try taking this and modifying it:

function RungeKuttaFourthOrderSingleVariable(f, a, h, n)
  x(1) = 0;
  y(1) = a;

  for j = 1 : (n/h)
     x(j+1) = (j*h);
     temp1 = h * f(x(j),y(j)); 
     temp2 = h * f(x(j) + 0.5*h, y(j) + 0.5*temp1);
     temp3 = h * f(x(j) + 0.5*h, y(j) + 0.5*temp2);
     temp4 = h * f(x(j) + h, y(j) + temp3);
     y(j+1) = y(j) + (1/6) * (temp1 + 2*temp2 + 2*temp3 + temp4);
  end
end

The f(x(j),y(j) represents the function handle for calculating the derivative given an x and a currently y value from where the derivative is being evaluated. The code is for evaluated one DE, but it can be easily modified to evaluate a system of DE's.

Also I should have asked, but forgot to: what is L_e a function of?

 

[RedAnsar]

Try taking this and modifying it:

function RungeKuttaFourthOrderSingleVariable(f, a, h, n)
  x(1) = 0;
  y(1) = a;

  for j = 1 : (n/h)
     x(j+1) = (j*h);
     temp1 = h * f(x(j),y(j)); 
     temp2 = h * f(x(j) + 0.5*h, y(j) + 0.5*temp1);
     temp3 = h * f(x(j) + 0.5*h, y(j) + 0.5*temp2);
     temp4 = h * f(x(j) + h, y(j) + temp3);
     y(j+1) = y(j) + (1/6) * (temp1 + 2*temp2 + 2*temp3 + temp4);
  end
end

The f(x(j),y(j) represents the function handle for calculating the derivative given an x and a currently y value from where the derivative is being evaluated. The code is for evaluated one DE, but it can be easily modified to evaluate a system of DE's.

Also I should have asked, but forgot to: what is L_e a function of?

Thanks! I'll try and work with that.
L_e is just the Langevin function, coth(x)-(1/x), with the variable ξ_e as the parameter: L_e = L(ξ_e) = coth(ξ_e) - 1/(ξ_e)

 

[RedAnsar]

Playing some more, I don't think the explicit function for RungeKutta is necessary. ode45 is a fourth-order RungeKutta method, after all.
I wrote some more code, this time including the second equation:
firstode.m

function Eprime = firstode(t,E)
Eprime=(((0.0016).^(-1))*((coth(E)-(1./E))/((E).^(-2)-(csch(E)).^(2)))*(((10)/E)*cos(625*t)-1));

secondode.m

function Fprime = secondode(t,F)
Fprime = (((0.0016).^(-1))-(1/(2*0.0016))*((1/(coth(E)-(E).^(-1)))-((E).^(-1)))*10*F*cos(625*t));

odesolver.m

E0=10;
sol_1=ode45(@firstode,[0 1],E0);
t = linspace(0,1,10);
E = deval(sol_1,t,1);
sol_2=ode45(@secondode,[0 1],E0);
F = deval(sol_2,t,1);

So, my thought process is:
- I get MATLAB to evaluate the first ODE. Now I have an array of numerical solutions in the interval given. So far so good.
- I need to somehow plug in this array of numerical values into the second ODE, as the solution to the first ODE, the variable E, appears in the second ODE.
- I want to somehow tell MATLAB: Hey, that E in the secondode M-file is the solution of the ODE you just solved. Use that.
However, I get an error when running into the last two lines of the the odesolver M-file
At first, I thought I could just rewrite
sol_2=ode45(@secondode,[0 1],E0);
into
sol_2=ode45(@secondode,[0 1],E0, E);
as if the values of E were some extra value I can just latch on, but that didn't work out.

Any ideas...?