I want to solve this system of differential equations. Is there any code that can solve this?
2 views (last 30 days)
Show older comments
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;
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1496466/image.jpeg)
0 Comments
Accepted Answer
Torsten
on 29 Sep 2023
Edited: Torsten
on 29 Sep 2023
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])
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
3 Comments
Torsten
on 29 Sep 2023
Edited: Torsten
on 30 Sep 2023
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)
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)])
More Answers (0)
See Also
Categories
Find more on Logical in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!