Problems solving cupled 2nd Order ODE with od45

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

[Jan]

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]

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]

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

Best wishes

Torsten.

[Erik Kostic]

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]

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.

[Erik Kostic]

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]

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]

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]

[t,y] = ode45(M,Interval,Conditions);

plot(y(:,1),y(:,3));

Best wishes

Torsten.

[Erik Kostic]

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

[Steven Lord]

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.

Answer 2

[Dariusz Skibicki]

Replace

V = odeToVectorField(ode1);

with

V = odeToVectorField(odes);

[Dariusz Skibicki]

And

Conditions = [0 0];

with

Conditions = [0 1 0 1];