Solving Coupled ODE's by ODE45

Dear all,

I have a coupled ODE with two functions. One of them is a gaussian beam. I do want to solve the other one. So far to understand how ode45 works I wrote down such a code which is given below.

 gamma = (1e+9)/3
omega = (1e+9)/3
[t,y] = ode45(@(t,y) arbit2(t,y,gamma,omega), [0 20e-9],[0 1]);
plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of coupled ode with ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2') 
          function dydt = arbit2(t,y, gamma, omega)
dydt = [(-gamma-omega)*y(1)+(omega)*y(2); (gamma+omega)*y(1)+(-omega)*y(2)];
          end

My question is how can I modify y(1) as a function. For example I want to write a function that gives y(1) as following;

 t = 0:1e-15:1e-12;
a=3
y = a*gaussmf(t,[0.5e-13 5e-13]);

I want to solve y(2) according to given y(1) Thanks in advance.

Answer 1

[Star Strider]

If I understand correctly what you want, you can simply modify ‘arbit’ (and calls to it) directly:

arbit2 = @(t,y, gamma, omega,a) [a*gaussmf(t,[0.5e-13 5e-13]); (gamma+omega)*y(1)+(-omega)*y(2)];

gamma = (1e+9)/3
omega = (1e+9)/3
a = 3;
tv = 0:1e-15:1e-12;
[t,y] = ode45(@(t,y) arbit2(t,y,gamma,omega,a), tv,[0 1]);
plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of coupled ode with ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')

It is not an interesting result.

[Ömer Yaman]

Actually, I try to find fluorescence time by introducing one of the coupled function as a laser pulse which has gaussian intensity profile. Let's call that function as y(1), with the relation of coupled ode I want to solve y(2).

As an input, I want to make y(1) gaussian as I stated before. I hope this explanation is clear.

[Star Strider]

Note that ‘y(1)’ is the first row of your ODE function. When I substituted ‘a*gaussmf(t,[0.5e-13 5e-13])’ for it and tweaked ‘arbit2’ to accept ‘a’ as an argument, your function works with ode45.

Unfortunately, it is not clear. I have no idea what you are doing with respect to finding ‘fluorescence time’, since that is not an area of my expertise.

EDIT — (12:48 UCT)

If you want ‘(a*gaussmf(t,[0.5e-13 5e-13])’ as an input to your system, this could work:

arbit2 = @(t,y, gamma, omega,a) [((-gamma-omega)*y(1)+(omega)*y(2)) + (a*gaussmf(t,[0.5e-13 5e-13])); (gamma+omega)*y(1)+(-omega)*y(2)];

It adds the Gaussian as a forcing function to ‘y(1)’.

[Ömer Yaman]

Actually I'm trying to solve that.

x'(t)=[(-a-b)*x(t)]+[b*y(t)] y'(t)=[(a+b)*x(t)]+[-b*y(t)]

x(t)= gaussian. I want to solve y(t)

[Star Strider]

I have no idea what

x(t)= gaussian

implies, or what ‘a’ and ‘b’ are.

Try this:

arbit2 = @(t,z,a,b) [(-a-b).*(exp(-(a-z(1)).^2/b)) + b.*z(2); (a+b).*(exp(-(a-z(1)).^2/b))+(-b-z(2))];

Answer 2

[Torsten]

function dydt = arbit2(t,y, gamma, omega)
a=3;
y_fix=a*gaussmf(t,[0.5e-13 5e-13]);
dydt = [(-gamma-omega)*y_fix+(omega)*y(2); (gamma+omega)*y_fix+(-omega)*y(2)];
          end

[Ömer Yaman]

if I set a=0, it gives still y(1) result. how?