Solving a differential equation

Hello,

I am working at a mathematical problem and I wanted to solve it numerically through MATLAB. The equation has the form:

y'(t) = g(t,y(t)) - a y(t)

and

g(t,y(t)) = b sinh[(g(t-dt, y(t-dt))] y(t) + c

with a,b,c in R.

I already tried using the ode and dde function in MATLAB without luck. Therefore I hope you can help me with the implementation in the code.

[Matt Tearle]

Please show what you've done and explain what the problem ("without luck") is.

Also, am I correct in interpreting the equation as being defined implicitly (g(t,y) is a function of g(t-dt,y(t-dt)))? I don't really understand how that would work.

[John D'Errico]

Matt - that is a variation of a DDE, a delayed differential equation.

[Andrew Newell]

I think Matt has a point. As formulated, it looks like an infinite recursion.

[Andrew Newell]

Thanks to whoever corrected the spelling of the title!

Answer 1

[Andrew Newell]

If your equations contain only linear combinations of the derivatives, you could formulate it as a mass matrix problem M(t,y)y' = f(t,y) . For the two equations you have provided in the comment,

dF = - R*dI/u_a;
dI = const2*cosh(insinh*F)*insinh*dF+vor*sinh(insinh*F)*dn*exp(-beta*dEC2);

you could set up the problem as follows:

y = [F; I];

M = [1 - R/u_a; const2*cosh(insinh*F)*insinh 1];

and f(t,y) would output

[0; vor*sinh(insinh*y(1))*dn*exp(-beta*dEC2)]

See the ode45 documentation for how to include a mass matrix.

(edited to include the relevant part of betlor5's comment)

Answer 2

[Andrew Newell]

I have been floundering a bit with this problem, but now I see what you need: a variable

z = [y; g]

with

y'(t) = (g-a)*y(t)

g'(t) = df(t)/dt * y(t)

g(0) = c

The function f is something like what you currently would call b sinh[(g(t,y)].

Note that this new system is a set of ordinary differential equations, so you could use the odexx functions.

EDIT: In terms of z,

z_1'(t) = (z_2-a)*z_1

z_2'(t) = df(z_1)/dt * z_1

z_2(0) = c

Answer 3

[betlor5]

The code is:

function d = ddeg(t,y,Z)
  A = 1000;
  dz = 6;
  N_T1 = 10^(-2);
  N_T2 = 10^(-2);
  dEC1 = 0.35; 
  dEC2 = 0.01; 
  tau_0 = 1e-15;
  tau_n = 5e-9; 
  B = 1; 
  V_A = 1.2;
  q = 1.602176487*10^(-19);
  R = 4e3;

    beta = 1/(8.617343*10^-5 * 298.15);

    insinh = beta*dz/2;
    vor = 2*q*A*dz/tau_0;

    d = [0;0];

    n_1 = N_T1;
    n_2 = N_T2 * exp(-(dEC1-dEC2)*beta);
    const = vor*n_1*exp(-beta*dEC1)+vor*n_2*exp(-beta*dEC2);
    I = (const + Z(1,1) * exp(-beta*dEC2))*sinh(insinh*(V_A-(R*Z(2,1))));
    F_OFF = asinh(I/const)/insinh;
    G = n_1*exp(-B/F_OFF)/tau_0;
    d(1) = G - y(1)/tau_n;
end

with using

sol = dde23('ddeg',[1e-15],[0,1.2e-6],[1e-15,1e-8]);

I get a completly diffrent answer than expected. I want to simulate it in an intervall from 1e-15 to 1e-8 with a minimal dt = 1e-15 between an interration. for y(1) the start value is 0 and for y(2)=I the start value is 1.2e-6.

I hope this helps solving the problem

Answer 4

[Andrew Newell]

If I understand your question, you have a delay differential equation with one variable y(t), and your equation has the form y'(t) = f( y(t), y(t-t1) ), so it should input a scalar for y(t) and another for y(t-1) (the variables y and Z). Instead, you input vectors of length 2 and return a 2-vector for y'(t).

EDIT: For clarity, it might help to rewrite your function as follows:

function d = ddeg(~,y,ylag)

The tilde is used in recent versions of MATLAB to indicate that the variable is not used. If your version of MATLAB objects, put the t back in.

Also, look through Matlab dde23 Tutorial M-files for examples of how to use dde23 properly.

EDIT 2: I'm not sure that your definition of the function g makes sense. If

g(t,y(t)) = b sinh[(g(t-dt, y(t-dt))] y(t) + c

then

g(t-dt,y(t-dt)) = b sinh[(g(t-2*dt, y(t-2*dt))] y(t-dt) + c

and so on forever. How could you evaluate g?

[betlor5]

If I use a scalar function for the dde solver, I will have no possibility to pass g(t,y(t)) and it's history to the next iteration. Therefore the problem still remains or you could ask how can I pass specific variables to the dde in order to use it in the next iteration step?

[Andrew Newell]

You don't need to pass on g. The purpose of ddeg is to calculate y'(t) and get the next y(t). The solver takes care of the details.

[betlor5]

How do I express I, if I am only using sclar functions. Right know I use:

I = (const + Z(1,1) * exp(-beta*dEC2))*sinh(insinh*(V_A-(R*Z(2,1))));

How do I express I(t-dt), which would be Z(2,1).

Answer 5

[betlor5]

It is a gerneral problem with recusive function that they will go one forever if you don't give them a start parameter.

The idear is that you say "g(t,y(t)) = b sinh[(g(t-dt, y(t-dt))] y(t) + c

then

g(t-dt,y(t-dt)) = b sinh[(g(t-2*dt, y(t-2*dt))] y(t-dt) + c

and so on" but if you at the time t= 0 you say

G(t=0,y(t=0)=0) = c

and in the code I defined c=1.2e-6. Inorder to get a better understending it would be rewritten as:

g = {g(t,y(t)) = b sinh[(g(t-dt, y(t-dt))] y(t) + c, t>0 and c, t =0

and the implementation of this parameter and it's recursive behaviour is the problem I have got in the mentioned code.

[betlor5]

We are only loking at an specific intervall.

Let's say t in [0,a].Therefore you would start at f(0) = c and have your value for

f(dt) = f(0) = c

In a mathematical way you would solve this equation through using a boundary condition. The same necessety is given for an differential equation. Therefore I thought using a second boundary condition in order to solve my two equations from above. Let's go back to my equation from above:

y'(t) = g(t,y(t)) - a y(t)

and

g(t,y(t)) = b sinh[(g(t-dt, y(t-dt))] y(t) + c

I thought it goes without saying that you use y(t=0) = y_0 and g(t=0,y_0) = g_0 as the boundary conditions. In the code that would be [0,1.2e-6]. At least that was my understending of how I have to implement the starting point for a de.

I hope this clarifies the question.

Answer 6

[betlor5]

Shouldn't g'(t) = f'(t) * y'(t) ?

So there is no way of passing g as a variable without using this solution?

The problem is that in the not simplified equation I would get a system of 8 differential equations. I allready tried solving that one, but it is way to complex in order to get an answer in an acceptable time.

[Andrew Newell]

It could be g'(t) = f'(t) * y'(t). My answer above is just a sketch; I can't tell exactly what the right form is without the original equations. Yes, I think that if you can calculate g'(t), you should make g a component of the y vector.

I don't see the problem with solving 8 differential equations. People solve systems of thousands using MATLAB. And I really don't see how passing g as a variable makes it any simpler.

[betlor5]

Here is an example of my code:

dF = - R*dI/u_a;

dI = const2*cosh(insinh*F)*insinh*dF+vor*sinh(insinh*F)*dn*exp(-beta*dEC2);

Therefore I have got the problem that I need dF in order to calculate dI and dI in order to calculate dF. How can you tell matlab to solve such an system of differential equation? Do I ignore the warning?