Problems solving cupled 2nd Order ODE with od45

2017 11 28
2 个回答
更新时间：2023 3 16
3 次查看（30 天）

0 个投票

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更早的评论

Jan
Jan 2017-11-28
This seems to be a homework problem. Then show us, what you have tried so far and ask a specific question.
Did you convert the equation of the 2nd order to a system of equations of 1st order already?
Erik Kostic
Erik Kostic 2017-11-29
Hello Jan,
yes, i have used the odeToVectorField method to convert both of the two eqations. And there is the problem, because i can't use from the new vector that was generated a Y(1) to solve via ode45.
Torsten
Torsten 2017-11-29
Why don't you just show what you have so far ?
Best wishes
Torsten.
Erik Kostic
Erik Kostic 2017-11-29
编辑：Torsten 2017-11-29
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);

请先登录，再进行评论。

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

回答（2 个）

Torsten
Torsten 2017-11-29

0 个投票

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
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);
Torsten
Torsten 2017-11-29
Delete the blank at the end of the definition of M:
Replace
...(mag3(3))^2)^3)) -C*y(3)]
with
...(mag3(3))^2)^3))-C*y(3)]
Best wishes
Torsten.
Erik Kostic
Erik Kostic 2017-11-29
编辑：Erik Kostic 2017-11-29
Hello, now it solves the ode perfectly.
Can i use
fplot(@(x)deval(Solution,x,1),[0 20])
to plot the behavoir of the pendulum in the plane? Again thank you very much
to be more precise, what i want to plot is the trajectory, from the starting position to the attraction point.
Torsten
Torsten 2017-11-29
[t,y] = ode45(M,Interval,Conditions);
plot(y(:,1),y(:,3));
Best wishes
Torsten.
Erik Kostic
Erik Kostic 2017-11-29
Hey Torsten, thank you very much you are a germ :D
Steven Lord
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.

请先登录，再进行评论。

Dariusz Skibicki
Dariusz Skibicki 2023-3-16

0 个投票

Replace
V = odeToVectorField(ode1);
with
V = odeToVectorField(odes);

类别

在 帮助中心 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

产品

标签

Community Treasure Hunt

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

Start Hunting!

Translated by