I want to solve this system of differential equations. Is there any code that can solve this?
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Akzeptierte Antwort
Torsten
am 29 Sep. 2023
Bearbeitet: Torsten
am 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 Kommentare
Torsten
am 29 Sep. 2023
Bearbeitet: Torsten
am 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)])
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Logical finden Sie in Help Center und File Exchange
Produkte
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!