how to solve a set of differential equation systems like this?

Dear friend,

How to solve a ordinary differential equation systems like above using MATLAB?

a=b=c=1

Your help would be highly appreciated.

[Torsten]

What's so special about the system that you think it cannot be solved by one of the standard MATLAB integrators (ODE45, ODE15S etc.) ?

[Sam Chak]

Hi @Daniel Niu,

You forgot to give the value of d on the right side of the ODEs. Can we assume that as well?

Also, at some time t, the denominator may cause singularity if , that is when and .

[John D'Errico]

@Sam Chak

Actually, the point where you identify a singularity is an interesting one, since in order for a singularity to exist, one would need to have

x^2 + (c*t-y)^2 == 0

And that implies two things MUST happen, x==0, AND y = c*t. Note that both differential equations have a corresponding term in the denominator. But then look carefully. In the numerator of both equations, they have a similar term on top. As such, that fraction is actually going to have a finite limit at that point. At the exact point of interest, the result is effectively a Nan, but everywhere else, it is well defined. Not really any different from the fraction x/x.

[Sam Chak]

@John D'Errico, Thank you for pointing out this and the explanation.

Answer 1

[Torsten]

ar = 0.01;
af = 0.001;
x0 = 3;
y0 = -3;
Lr = 0;
Lf = 0;
z0 = [x0; y0; Lr; Lf];
[T,Z] = ode45(@(t,z)fun(t,z,ar,af),[0 11.5],z0);
plot(Z(:,1),Z(:,2),Z(:,3),zeros(size(T)))
ylim([-3 0.5]) 

function dz = fun(t,z,ar,af)
  x = z(1);
  y = z(2);
  Lr = z(3);
  Lf = z(4);
  sr = exp(-ar*Lr);
  sf = exp(-af*Lf);
  dz = zeros(4,1);
  dz(1) = sf* (Lr-x)/sqrt((Lr-x)^2+(0-y)^2);
  dz(2) = sf* (0-y)/sqrt((Lr-x)^2+(0-y)^2);
  dz(3) = sr;
  dz(4) = norm([dz(1) dz(2)]);
end

[Daniel Niu]

Dear Torsten,

Sorry to bother you again. I want to know if the rabbit is not moving horizontal but from (0,0) to (-100,100)

how to modify the code?

Thank you!

[Torsten]

Uninteresting case since both fox and rabbit will run on the same straight line.

ar = 0.01;
af = 0.001;
x0r = 0;
y0r = 0;
x0f = 3;
y0f = -3;
Lr0 = 0;
Lf0 = 0;
z0 = [x0r; y0r; x0f; y0f; Lr0; Lf0];
[T,Z] = ode45(@(t,z)fun(t,z,ar,af),[0 11.5],z0);
plot(Z(:,1),Z(:,2),Z(:,3),Z(:,4))

function dz = fun(t,z,ar,af)
  xr = z(1);
  yr = z(2);
  xf = z(3);
  yf = z(4);
  Lr = z(5);
  Lf = z(6);
  sr = exp(-ar*Lr);
  sf = exp(-af*Lf);
  dz = zeros(6,1);
  dz(1) = sr* (-1/sqrt(2));
  dz(2) = sr* (1/sqrt(2));
  dz(3) = sf* (xr-xf)/sqrt((xr-xf)^2+(yr-yf)^2);
  dz(4) = sf* (yr-yf)/sqrt((xr-xf)^2+(yr-yf)^2);
  dz(5) = norm([dz(1) dz(2)]);
  dz(6) = norm([dz(3) dz(4)]);
end

Answer 2

[Torsten]

Hint:

The length L(t) of a curve in parametric form (x(t),y(t)) where x(t) and y(t) are given by differential equations

dx/dt = f(x,y)
dy/dt = g(x,y)

can be computed by additionally solving a differential equation for L:

dL/dt = sqrt(f(x,y).^2+g(x,y).^2)

with initial condition

L(0) = 0.

[Daniel Niu]

Dear Torsten,

Thank you for your answers. I need more help from you.

The key problem here is the speed sf is variable, and sf=exp(-a*df),where df is the distance it have travelled.

[Torsten]

Shouldn't there be an analogous equation for r1 as there is for z1 and z2 ?

dr1/dt = exp(-a*r1), r1(0) = 0

?

syms r1(t) a
sol = dsolve(diff(r1,t)-exp(-a*r1)==0,r1(0)==0)

sol = 

a = 0.001;
t = 0:1:10000;
plot(t,log(a*t+1)/a)

[Daniel Niu]

Dear Torsten,

the problem is sf is variable value and I want solve the ODE for z1 and z2.

I don't need the value of r1.

thank you!

[Torsten]

r1 is part of your ODE system:

dz1/dt = s_f* (r1-z1)/sqrt((r1-z1)^2+(r2-z2)^2)

dz2/dt = s_f* (r2-z2)/sqrt((r1-z1)^2+(r2-z2)^2)

And I think the speed of the rabbit will also decrease with time...

[Sam Chak]

Hi @Daniel Niu

Your original "Fox-chasing-Rabbit" example assumes that the Rabbit is moving at a constant velocity . That's why the solution for the Rabbit's position is .

If the Rabbit's velocity is time-varying, or is reactively dependent on the position of the Fox {, }, then naturally there exists , and I think that is what @Torsten tried to tell you.