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

조회 수: 1 (최근 30일)
Student 2023년 9월 29일
댓글: Student 2023년 10월 3일
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 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이전 댓글 1개 표시이전 댓글 1개 숨기기
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 2023년 10월 3일
this code only works at t<0.4 can you solve this problem?(I want for t>10)

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

### 카테고리

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