Coupled system of ODE's using ode23

Dear reader,

I am currently trying to solve a system of ODE's using matlab. Sadly I seem to have a hard time understanding the function ode23 and the way one should use it for this kind of problem. The ODE's I am trying to solve for are:

Now according to the ode23 discription it needs a function input and so I made a .m file called ''DEdef.m'' and contains the following lines of code:

function [] = DEdef(I, D)
%IV,I,IVsol - Independent variables
%DV,D,DVsol - Dependent variables

%dependent veriables are N2 and q
N2=D(1);
q=D(2);        
Ddv_Div=[W_p*N-sigma*c*N2*q/V-N2/t_2;            
         q*(sigma*c*N2-1/t_2)]                   
end

I am then trying to solve in using a different .m file with the following lines of code:

clc
clear all

%constants
c=300000*10^-7;          %cm/ns
t_res=20;                %ns
t_2=230000;              %ns
sigma=4*10^-19;          %cm2
V=((2/10)^2*pi)*10;      %volume cm3
N=10^19;                 %cm-3
W_p=1;                   %no idea.....
N2_0=0.1;                %starting value
q0=0.1;                  %starting value


IC=[N2_0 q0];           %Initial values
t=[0:1:10*t_2];         %time domain   


[IVsol,DVsol] = ode23('DEdef' , domain, IC);

sadly this results in an error stating that there are ''too many output arguments''. Is there anyone with more experience using the ODE functions and would like to explain to me what I am doing wrong and what I might wanna do differently in order to get the script running?

Thanks in advance :)!

Answer 1

[Steven Lord]

There are a couple problems and at least one area for improvement for your code. I'll be referring to several sections of the documentation page for ode23 so I'll leave the link here for reference.

The two main problems for your attempt lie with your ODE function that evaluates the right-hand side of your system.

function [] = DEdef(I, D)
%IV,I,IVsol - Independent variables
%DV,D,DVsol - Dependent variables

%dependent veriables are N2 and q
N2=D(1);
q=D(2);        
Ddv_Div=[W_p*N-sigma*c*N2*q/V-N2/t_2;            
         q*(sigma*c*N2-1/t_2)]                   
end

The first problem is that your function does not satisfy the requirement given in the section of the documentation that talks about the odefun input: "The function dydt = odefun(t,y), for a scalar t and a column vector y, must return a column vector dydt of data type single or double that corresponds to f(t,y). odefun must accept both input arguments, t and y, even if one of the arguments is not used in the function."

Your function returns no output arguments. Modify your code so it returns Ddv_Div.

The second problem is that your function uses a lot of undefined variables. Yes, I know you defined them in your script file from which you called ode23 but that does not make them visible to your function. Either define them inside your DEdef function (my preferred approach) or pass them into DEdef as additional parameters.

In your script file where you call ode23:

[IVsol,DVsol] = ode23('DEdef' , domain, IC);

you're using a variable you haven't defined (domain, which in context should probably be t.) As one area for improvement, you're using a very old syntax for the ODE solvers. While that syntax, where you pass a char vector containing the name of your ODE function, is still accepted for backwards compatibility I strongly recommend you specify a function handle as the first input to the solver as the examples on the documentation page (especially the "Solve Nonstiff Equation" example) do.

[Jesse Verstappen]

Dear Steven,

First and foremost thanks for the quick and above all elaborate response to my question.

I did indeed notice that my function format does not correspond to the format used by MATLAB, thanks for pointing that out! I also included the parameters in the function format as you stated that this might be important (I only need to change the W_p anyways and I can always add this additional parameter thing after getting the script working). I now have the DEdef.m as follows:

function dydt=DEdef(t,y)

%Constants:
c=300000;                 %m/s
t_res=20*10^-9;          %s
t_2=230*10^-6;           %s
sigma=4*10^-23;          %m2
V=pi*4*10^-7      ;      %volume m3
N=10^25;                 %m-3

N2c=1/(c*sigma*t_res);
W_pc=N2c/(N*t_2);

W_p=W_pc;

%ODE definition
dydt=[0 ; 0];
dydt(1)=W_p*N-sigma*c*y(1)*y(2)/V-y(1)/t_2;
dydt(2)=y(2)*(sigma*c*y(1)-1/t_2);                  
end

Please note that I did rewrite most of the constants to get a phyicsal result. I then defined the time vector as:

number_of_points=100;
t_step=(0.5*t_2-t_start)/(number_of_points-1);
t_domain=[t_start:t_step:0.5*t_2]; 

Where t_start is zero and the step function simply creastes an evenly spaced time interval.

The ode23 function now looks like the following line of code:

[t,y] = ode23(@DEdef , t_domain, IC);

However even for 100 steps this seems to get stuck (it never gives a final result no matter how long I wait). Now just to be clear I do know what kind of solution this should give and I will add those in the attachment, as this might give you some more insight in what is going on. I do expect the computation to be rather quick due to the low ammount of steps i am forcing it to do (yes I know the poor resolution might be giving an instable result, but for now I am trying to get any result at all, suck that I can up scale the model lateron to 10k time steps).

As the script does run, but never spits out any kind of result I am guessing I am still missing something rather important. Could you maybe help me out once more?

Thanks in advance!

With kind regards,

Jesse

[Steven Lord]

Let's look at the constants you're using.

c=300000;                 %m/s
t_res=20*10^-9;          %s
t_2=230*10^-6;           %s
sigma=4*10^-23;          %m2
V=pi*4*10^-7      ;      %volume m3
N=10^25;                 %m-3
N2c=1/(c*sigma*t_res)

N2c = 4.1667e+24

W_pc=N2c/(N*t_2)

W_pc = 1.8116e+03

W_p=W_pc;

Both N and N2c are quite large, while t_res, V, and (especialy) sigma are very small. Let's look at what your equations are when we use those constants.

syms y1(t) y2(t)
y = [y1(t), y2(t)];
dydt=sym([0 ; 0]);
dydt(1)=W_p*N-sigma*c*y(1)*y(2)/V-y(1)/t_2;
dydt(2)=y(2)*(sigma*c*y(1)-1/t_2);
vpa(dydt, 5)

ans = 

Unless the components of the solution are very large, the first term in dydt(1) dominates and y1 takes off like a rocket. That will counteract the small coefficient in dydt(2) and cause it to rise rapidly as well, which could satisfy the "Unless the components of the solution are very large" from my previous sentence.

When I run the attached example, which is mostly your code, I see a solution with one component that rockets up to around 3.25e23 then crashes down to around 1e21 by x = 1.15e-4. The other component of the solution is "only" 2e18 at that point.

You might want to double-check your dimensionals and you may also want to try a stiffer solver.

[Jesse Verstappen]

Dear Steven,

Thanks once more for the reply. The dimensions are correct. I also expect the solution of the first ODE to jump up and down, however it should get stable after a some point has past. the second ODE should jump up and down as well, as it is zero below the critical value of the first ODE and nonzero above.

The large values are due to the physics behind the ODE. The set of ODE are used to simulate the startup phase of a laser and give the inverse population density and photon count respectively. Ofcourse only a small volume (V) of gain medium results in large amounts of e- being excited (N2) (and falling back after some time, ''creating'' a photon (q)).

I do have 1 more question with regards to MATLAB and solving ODEs. As after lots of extra work on the model and some help from another physics student I actually got the script running:

%Constants:
c=299792458;                 %m/s
t_res=20*10^-9;          %s
t_2=230*10^-6;           %s
sigma=4*10^-23;          %m2
V=pi*4*10^-7      ;      %volume m3
N=10^25;                 %m-3

N2c=1/(c*sigma*t_res);
W_pc=N2c/(N*t_2);

%W_p=W_pc;
W_p=1.2*W_pc;
q0=(W_p-W_pc)*N*V*t_res;


A=W_p*N;
B=sigma*c/V;
C=t_2;
D=sigma*c;
E=1/t_res;

t_domain=[0:t_2/1000:10*t_2];
%%
dydt=@(t,y)[A-B*y(1)*y(2)-y(1)/C;y(2)*(D*y(1)-E)]
[t,y]=ode23(dydt,t_domain,[0.1 0.1], odeset('OutputFcn', @odeplot))
figure(2)
subplot(1,2,1)
plot(t,y(:,1))
subplot(1,2,2)
plot(t,y(:,2))

As you can see I only use the dydt@ no functions.m files are used and it only contains 2 lines of code that do the hard work. This leaves me wondering why the method explained by the ode23-page actually seems to be unable to sovle it and this simple code does and why this type of solution is not included in the ode23 description?

Thanks in advance!