필터 지우기
필터 지우기

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

조회 수: 2 (최근 30일)
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;

채택된 답변

Torsten
Torsten 2023년 9월 29일
편집: Torsten 2023년 9월 29일
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
  댓글 수: 3
Torsten
Torsten 2023년 9월 29일
편집: Torsten 2023년 9월 30일
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
Student 2023년 10월 3일
this code only works at t<0.4 can you solve this problem?(I want for t>10)

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Logical에 대해 자세히 알아보기

제품

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by