I want to solve this system of differential equations. Is there any code that can solve this?

I am trying to solve this system of differential equations using numerical analysis. However, I couldn't find any Matlab code to solve this.

Is there a code that can numerically analyze this type of equation?

I = 0.000122144; %kgm^2
m = 0.19744; %kg
g = 9.80665; %m/s^2
R = 0.035; %m
s = 0.090757; %rad
r = -0.011515; %m
M = 0.029244; %kg
k = 15.36;
A = 0.011065331; %m^2
syms x(t) y(t)
ode1 = I * diff(y, 2) == -m * sqrt(g^2 + R^2 * (diff(x, 2))^2 - 2 * g * R * diff(x, 2) * sin(s)) * r * sin(s + y - atan(R * diff(x, 2) * cos(s) / (g - R * diff(x, 2) * sin(s)))) - k * R * A * (diff(y) - diff(x));
ode2 = 2 * M * R * diff(x, 2) == (M + m) * g * sin(s) - m * r * (sin(y) * (diff(y)^2) - cos(x) * (diff(x)^2));
condition1 = x(0) == 0;
condition2 = y(0) == 0;
condition3 = diff(x) == 0;
condition4 = diff(y) == 0;

Answer 1

[Torsten]

Note that there is a discrepancy between your code and your graphics for the second ode. I took the version from the graphics ( which worked better for ode15i :-) )

y0 = [0 0 0 0];
yp0 = [0 0 0 0];
t0 = 0;
tspan = [t0,0.25];
[y0_new,yp0_new] = decic(@f,t0,y0,[1 1 1 1],yp0,[0 0 0 0])

y0_new = 4×1

0 0 0 0

yp0_new = 4×1

0 0 -172.0836 289.5405

sol = ode15i(@f,tspan,y0_new,yp0_new,odeset('RelTol',1e-8,'AbsTol',1e-8));
figure(1)
plot(sol.x,[sol.y(1,:);sol.y(2,:)])

[y,yp] = deval(sol,sol.x);
for i = 1:numel(sol.x)
  res(i,:) = f(sol.x(i),y(:,i),yp(:,i));
end
figure(2)
plot(sol.x,[res(:,3),res(:,4)])

function res = f(t,y,yp)
  I = 0.000122144; %kgm^2
  m = 0.19744; %kg
  g = 9.80665; %m/s^2
  R = 0.035; %m
  s = 0.090757; %rad
  r = -0.011515; %m
  M = 0.029244; %kg
  k = 15.36;
  A = 0.011065331; %m^2
  res = zeros(4,1);
  res(1) = yp(1) - y(3);
  res(2) = yp(2) - y(4);
  res(3) = I*yp(3) - (-m * sqrt(g^2 + R^2 * yp(4)^2 - 2 * g * R * yp(4) * sin(s)) * r * sin(s + y(1) - atan(R * yp(4) * cos(s) / (g - R * yp(4) * sin(s)))) - k * R * A * (y(3) - y(4)));
  res(4) = 2 * M * R * yp(4) -( (M + m) * g * sin(s) - m * r * (sin(y(1)) * (y(3)^2) - cos(y(1)) * yp(3)));
end

[Student]

Many thanks to you! I really owe you a lot.

[Torsten]

Better also look at the residuals to check whether the equations are really satisfied for the given solution. I changed the code above accordingly.

Maybe the code performs more reliable with analytical Jacobians:

I = 0.000122144; %kgm^2
m = 0.19744; %kg
g = 9.80665; %m/s^2
R = 0.035; %m
s = 0.090757; %rad
r = -0.011515; %m
M = 0.029244; %kg
k = 15.36;
A = 0.011065331; %m^2

syms t
y = sym('y',[4 1]);
yp = sym('yp',[4 1]);
f(1,1) = yp(1) - y(3);
f(2,1) = yp(2) - y(4);
f(3,1) = I*yp(3) - (-m * sqrt(g^2 + R^2 * yp(4)^2 - 2 * g * R * yp(4) * sin(s)) * r * sin(s + y(1) - atan(R * yp(4) * cos(s) / (g - R * yp(4) * sin(s)))) - k * R * A * (y(3) - y(4)));
f(4,1) = 2 * M * R * yp(4) -( (M + m) * g * sin(s) - m * r * (sin(y(1)) * (y(3)^2) - cos(y(1)) * yp(3)));
Jac_y = jacobian(f,y);
Jac_yp = jacobian(f,yp);
f = matlabFunction(f,'Vars',{t,y,yp});
Jac_y = matlabFunction(Jac_y,'Vars',{t,y,yp});
Jac_yp = matlabFunction(Jac_yp,'Vars',{t,y,yp});
Jac = @(t,y,yp) deal(Jac_y(t,y,yp),Jac_yp(t,y,yp));

y0 = [0 0 0 0].';
yp0 = [0 0 0 0].';
t0 = 0;
tspan = [t0,0.25];
options = odeset('RelTol',1e-8,'AbsTol',1e-8,'Jacobian',Jac);
[y0_new,yp0_new] = decic(f,t0,y0,[1 1 1 1].',yp0,[0 0 0 0].',options)

y0_new = 4×1

0 0 0 0

yp0_new = 4×1

0 0 -172.0836 289.5405

sol = ode15i(f,tspan,y0_new,yp0_new,options);
figure(1)
plot(sol.x,[sol.y(1,:);sol.y(2,:)])

[y,yp] = deval(sol,sol.x);
for i = 1:numel(sol.x)
  res(i,:) = f(sol.x(i),y(:,i),yp(:,i));
end
figure(2)
plot(sol.x,[res(:,3),res(:,4)])

[Student]

this code only works at t<0.4 can you solve this problem?(I want for t>10)