How to solve and plot system of nonlinear differential equations?

Hello,

I am having troubles solving a system of second order nonlinear equations with boundary conditions.

The equations are:

u0=4*pi*10e-7;
R=0.8;
V=10;
d=0.015;
xp1=0.002;
xp2=d-xp1;
p2=1000000;
p1=1000000;
r=0.01;
W=0.02;
g=0.001;
N=200;
x0=0.025;
k=300;
M=0.3;
B=40;

syms t x(t) i(t)
v(t) = 10*heaviside(0.3-t);
fplot(v(t), [0 0.6])

L0=u0*pi*r^2*N^2;
L=L0/(x+(r*g/(2*W)));

Eqn1:

The equation 2 has boundries, I can write in Mathmatica as,

-M*Dt[x[t], {t, 2}] - B*Dt[x[t], t] + 
  If[x[t] > xp2, -p2*(x[t] - xp2), 0] +
  If[x[t] < xp1, -p1*(x[t] - xp1), 0] - k*(x[t] - x0) + 
  0.5*(i[t]*i[t]*D[L[x], x[t]]) = 0, 
x[0] = xp2 - 0.00000001, x'[0] = 0., i[0] = 0.0

but I dont know how to write the second equation in matlab and solving and plotting for i(t) and x(t).

[Torsten]

I don't see any correspondence between the equation you posted and the equation you set up in Mathematica.

[berk26092]

Yes, they are 2 different equations. I am trying to solve a system of nonlinear equations.

I also wrote the first one in mathmatica as,

R*i[t] + L[x]*Dt[i[t], t] + i[t]*D[L[x], x[t]]*Dt[x[t], t] = V0[t]

[Torsten]

And where is the equation for L ?

[berk26092]

Okay, my bad. Now, I included in all detail.

[Torsten]

p1, p2, xp1, xp2, k ?

[berk26092]

Sorry... When I coded in mathmatica I used different names for them and in matlab I used different names. I think now the question should be all clear.

[Torsten]

When I coded in mathmatica I used different names for them and in matlab I used different names.

Why ? Don't make life harder than it needs to be.

Answer 1

[Torsten]

d=0.015;
xp1=0.002;
xp2=d-xp1;
tspan = [0 0.17];
y0 = [xp2 - 0.00000001;0;0];
options = odeset('RelTol',1e-12,'AbsTol',1e-12);
[T,Y] = ode15s(@fun,tspan,y0,options);
figure(1)
plot(T,Y(:,1))

figure(2)
plot(T,Y(:,2))

figure(3)
plot(T,Y(:,3))

function dy = fun(t,y)
  u0=4*pi*10e-7;
  r=0.01;
  N=200;
  g=0.001;
  W=0.02;
  M=0.3;
  B=40;
  R = 0.8;
  d=0.015;
  xp1=0.002;
  xp2=d-xp1;
  p2=1000000;
  p1=1000000;
  k=300;
  x0=0.025;
  L0=u0*pi*r^2*N^2;
  x = y(1);
  xp = y(2);
  i = y(3);
  L=L0/(x+(r*g/(2*W)));
  dLdx = -L0/(x+(r*g/(2*W)))^2;
  ToAdd = 0.0;
  if x < xp1
    ToAdd = -p1*(x - xp1);
  elseif x > xp2
    ToAdd = -p2*(x - xp2);
  end
  V0 = 0.0;
  if t <= 0.3
    V0 = 10;
  end
  dy = zeros(3,1);
  dy(1) = y(2);
  dy(2) = (-B*y(2) + ToAdd - k*(x-x0) + 0.5*i^2*dLdx)/M;
  dy(3) = (V0 - i*dLdx*xp - R*i)/L;
end

[berk26092]

Thank you!