Problems solving cupled 2nd Order ODE with od45

Question

Erik Kostic 2017-11-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/369778-problems-solving-cupled-2nd-order-ode-with-od45

评论： Dariusz Skibicki 2023-3-16

Hello.

I am given the task of simulating the two-dimensional motion of a magnetic pendulum in the x-y-plane. The problem comes down in solving this system of cupled 2nd order ordinary differential equation:

x'' + R*x' + sum_{i=1}^3 (m_i-x)/(sqrt((m1_i-x)^2 + (m2_i-y)^2 + d^2))^3 + G*x == 0

y'' + R*y' + sum_{i=1}^3 (m_i-y)/(sqrt((m1_i-x)^2 + (m2_i-y)^2 + d^2))^3 + G*y == 0

Those eqations discribe the motion in the plane. I know i can use the method "ode45" to solve such a problem, given some initial values.

I have tried it a few times, but didn't came to a solution.

I hope someone can help me. (x',y') = 0 no initial velocity and position (x,y) could be anywhere.

GREETINGS

4 个评论
显示 2更早的评论隐藏 2更早的评论

Torsten 2017-11-29

Why don't you just show what you have so far ?

Best wishes

Torsten.

Erik Kostic 2017-11-29

编辑：Torsten 2017-11-29

在 MATLAB Online 中打开

Hello Torsten

clear all, clc;
%%Constants
R = 0.2;
C = 0.3;
d = 0.5;
a = 1;
%%Position of magnets with input a,d > 0
mag1 = [ a/2, -sqrt(3)*a, -d];
mag2 = [-a/2, -sqrt(3)*a, -d];
mag3 = [   0,  sqrt(3)*a, -d];
%%Position of mass
pmp = [-10, -15,  0];
%%Velocity of mass
pmv = [ 0,  0,  0];
%%Acceleration of mass
pma = [ 0,  0,  0];
%%Matrix of trajectories
PMPos = zeros(3,1);
PMPos(:,1) = pmp;
%%ODE Solving
syms x(t) y(t)
ode1 = diff(x,t,2) + R*diff(x,t,1) - ( (mag1(1)-x)/(sqrt((mag1(1)-x)^2+(mag1(2)-y)^2+(mag1(3))^2)^3) + ...
                                       (mag2(1)-x)/(sqrt((mag2(1)-x)^2+(mag2(2)-y)^2+(mag2(3))^2)^3) + ...
                                       (mag3(1)-x)/(sqrt((mag3(1)-x)^2+(mag3(2)-y)^2+(mag3(3))^2)^3) ) +C*x == 0;
ode2 = diff(y,t,2) + R*diff(y,t,1) - ( (mag1(2)-y)/(sqrt((mag1(1)-x)^2+(mag1(2)-y)^2+(mag1(3))^2)^3) + ...
                                       (mag2(2)-y)/(sqrt((mag2(1)-x)^2+(mag2(2)-y)^2+(mag2(3))^2)^3) + ...
                                       (mag3(2)-y)/(sqrt((mag3(1)-x)^2+(mag3(2)-y)^2+(mag3(3))^2)^3) ) +C*y == 0;
odes = [ode1; ode2];
V = odeToVectorField(ode1);
M = matlabFunction(V,'vars', {'t','Y'});
Interval = [0 20];
Conditions = [0 0];
Solution = ode45(M,Interval,Conditions);

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2017-11-29

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/369778-problems-solving-cupled-2nd-order-ode-with-od45#answer_293766

在 MATLAB Online 中打开

M=@(t,y)[y(2);-R*y(2)+((mag1(1)-y(1))/(sqrt((mag1(1)-y(1))^2+(mag1(2)-y(3))^2+(mag1(3))^2)^3)+(mag2(1)-y(1))/(sqrt((mag2(1)-y(1))^2+(mag2(2)-y(3))^2+(mag2(3))^2)^3)+(mag3(1)-y(1))/(sqrt((mag3(1)-y(1))^2+(mag3(2)-y(3))^2+(mag3(3))^2)^3) )-C*y(1);y(4);-R*y(4)+((mag1(2)-y(3))/(sqrt((mag1(1)-y(1))^2+(mag1(2)-y(3))^2+(mag1(3))^2)^3) +(mag2(2)-y(3))/(sqrt((mag2(1)-y(1))^2+(mag2(2)-y(3))^2+(mag2(3))^2)^3) +(mag3(2)-y(3))/(sqrt((mag3(1)-y(1))^2+(mag3(2)-y(3))^2+(mag3(3))^2)^3) ) -C*y(3)];
Interval=[0 20];
Conditions = [x; dx/dt; y ; dy/dt] at t=0 ??
Solution = ode45(M,Interval,Conditions);

Best wishes

Torsten.

6 个评论
显示 4更早的评论隐藏 4更早的评论

Erik Kostic 2017-11-29

Hello Torsten, thank you for your quick answer. So am i guessing right, that the code reduces to the following?

clear all, clc; %% Constants R = 0.2; C = 0.3; d = 0.5; a = 1;

%% Position of magnets with input a,d > 0 mag1 = [ a/2, -sqrt(3)*a, -d]; mag2 = [-a/2, -sqrt(3)*a, -d]; mag3 = [ 0, sqrt(3)*a, -d];

%% Position of mass pmp = [-10, -15, 0];

%% Velocity of mass pmv = [ 0, 0, 0];

%% Acceleration of mass pma = [ 0, 0, 0];

%% Matrix of trajectories PMPos = zeros(3,1); PMPos(:,1) = pmp;

%% ODE Solving

M=@(t,y)[y(2);-R*y(2)+((mag1(1)-y(1))/(sqrt((mag1(1)-y(1))^2+(mag1(2)-y(3))^2+(mag1(3))^2)^3)+(mag2(1)-y(1))/(sqrt((mag2(1)-y(1))^2+(mag2(2)-y(3))^2+(mag2(3))^2)^3)+(mag3(1)-y(1))/(sqrt((mag3(1)-y(1))^2+(mag3(2)-y(3))^2+(mag3(3))^2)^3) )-C*y(1);y(4);-R*y(4)+((mag1(2)-y(3))/(sqrt((mag1(1)-y(1))^2+(mag1(2)-y(3))^2+(mag1(3))^2)^3) +(mag2(2)-y(3))/(sqrt((mag2(1)-y(1))^2+(mag2(2)-y(3))^2+(mag2(3))^2)^3) +(mag3(2)-y(3))/(sqrt((mag3(1)-y(1))^2+(mag3(2)-y(3))^2+(mag3(3))^2)^3) ) -C*y(3)]; Interval=[0 20]; Conditions = [-1; 0; -1; 0]; Solution = ode45(M,Interval,Conditions);

Or am i missing something, because after compiling it says:

Error using vertcat Dimensions of matrices being concatenated are not consistent.

Error in @(t,y)[y(2);-R*y(2)+((mag1(1)-y(1))/(sqrt((mag1(1)-y(1))^2+(mag1(2)-y(3))^2+(mag1(3))^2)^3)+(mag2(1)-y(1))/(sqrt((mag2(1)-y(1))^2+(mag2(2)-y(3))^2+(mag2(3))^2)^3)+(mag3(1)-y(1))/(sqrt((mag3(1)-y(1))^2+(mag3(2)-y(3))^2+(mag3(3))^2)^3))-C*y(1);y(4);-R*y(4)+((mag1(2)-y(3))/(sqrt((mag1(1)-y(1))^2+(mag1(2)-y(3))^2+(mag1(3))^2)^3)+(mag2(2)-y(3))/(sqrt((mag2(1)-y(1))^2+(mag2(2)-y(3))^2+(mag2(3))^2)^3)+(mag3(2)-y(3))/(sqrt((mag3(1)-y(1))^2+(mag3(2)-y(3))^2+(mag3(3))^2)^3)),-C*y(3)]

Error in odearguments (line 87) f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 113) [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in sim (line 31) Solution = ode45(M,Interval,Conditions);

Erik Kostic 2017-11-29

Hey Torsten, thank you very much you are a germ :D

Steven Lord 2017-11-29

Consider specifying the 'OutputFcn' option in your ode45 call as part of the options structure created by the odeset function. There are a couple of output functions included with MATLAB (the description of the OutputFcn option on that documentation page lists them) and I suspect one of odeplot, odephas2, or odephas3 will be of use to you.

请先登录，再进行评论。