Help solving a system of differential equations

Question

Tejas Adsul 2018-5-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/400867-help-solving-a-system-of-differential-equations

评论： Tejas Adsul 2018-5-16

I have the following system of differential equations, and I am not able to understand the best way to go about them. I tried using a couple of functions like dsolve, ode45, etc., but most of them give errors that I am not able to understand.

x(1)=0.3; x(3)=0.9; y(1)=0.4; y(3)=0.8;
A=10; p=10;
syms X1 X2;
num = ((x(3)-X1)-(X1-x(1)))*X1 + ((y(3)-X2)-(X2-y(1)))*X2;
den = sqrt((x(3)-X1)-(X1-x(1))^2 + (y(3)-X2)-(X2-y(1))^2);
Em = -A*(X1-x(1))*(x(3)-X1) - A*(X2-y(1))*(y(3)-X2) + p*num/den;
dXdt = [-diff(Em,X1); -diff(Em,X2)];

I would like to get X1 and X2 as a function of time. If I define syms X1(t) and X2(t), I get the error 'All arguments, except for the first one, must not be symbolic functions.' The above code, when run, gives me expressions in terms of X1 and X2. I need solutions of X1 and X2 in terms of t, where dX1dt = -diff(Em,X1), dX2dt = -diff(Em,X2).

Any help is appreciated. Thank you!

2 个评论
显示无隐藏无

Stephan 2018-5-15

Hi,

is it correct, that x and y are depending on t --> so: x(t) and y(t)?

Best regards

Stephan

Tejas Adsul 2018-5-15

编辑：Tejas Adsul 2018-5-15

The x and y whose values have been mentioned are constants. The variables X1 and X2 do depend on time. However, what I can do is ignore this time dependence, find the diff(Em,X1) and diff(Em,X2) expressions, then consider the time dependence by writing out the ode diff(X1,t) = -diff(Em,X1) and diff(X2,t) = -diff(Em,X2). I was able to do the first part, and am stuck with the second part.

请先登录，再进行评论。

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

Answer 1

Stephan 2018-5-15

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/400867-help-solving-a-system-of-differential-equations#answer_320342

编辑：Stephan 2018-5-15

在 MATLAB Online 中打开

Hi,

does this work for your purpose?

x1=0.3;
x3=0.9; 
y1=0.4;
y3=0.8;
A=10; 
p=10;
syms X_1 X_2 X1(t) X2(t);
num = ((x3-X_1)-(X_1-x1))*X_1 + ((y3-X_2)-(X_2-y1))*X_2;
den = sqrt((x3-X_1)-(X_1-x1)^2 + (y3-X_2)-(X_2-y1)^2);
Em = -A*(X_1-x1)*(x3-X_1) - A*(X_2-y1)*(y3-X_2) + p*num/den;
dEm_dX_1 = diff(Em,X_1);
dEm_dX_2 = diff(Em,X_2);
dEm_dX_1 = (subs(dEm_dX_1, [X_1 X_2], [X1 X2]));
dEm_dX_2 = (subs(dEm_dX_2, [X_1 X_2], [X1 X2]));
ode1 = diff(X1,t) == -dEm_dX_1;
ode2 = diff(X2,t) == -dEm_dX_2;
ode = matlabFunction([ode1; ode2])

ode is a function handle:

ode =
    function_handle with value:
      @(t)[diff(X1(t),t)==X1(t).*-2.0e1+(X1(t).*4.0e1-1.2e1).*1.0./sqrt(-X1(t)-X2(t)-(X2(t)-2.0./5.0).^2-(X1(t)-3.0./1.0e1).^2+1.7e1./1.0e1)+(X1(t).*2.0+2.0./5.0).*(X1(t).*(X1(t).*2.0-6.0./5.0).*1.0e1+X2(t).*(X2(t).*2.0-6.0./5.0).*1.0e1).*1.0./(-X1(t)-X2(t)-(X2(t)-2.0./5.0).^2-(X1(t)-3.0./1.0e1).^2+1.7e1./1.0e1).^(3.0./2.0).*(1.0./2.0)+1.2e1;diff(X2(t),t)==X2(t).*-2.0e1+(X2(t).*4.0e1-1.2e1).*1.0./sqrt(-X1(t)-X2(t)-(X2(t)-2.0./5.0).^2-(X1(t)-3.0./1.0e1).^2+1.7e1./1.0e1)+(X2(t).*2.0+1.0./5.0).*(X1(t).*(X1(t).*2.0-6.0./5.0).*1.0e1+X2(t).*(X2(t).*2.0-6.0./5.0).*1.0e1).*1.0./(-X1(t)-X2(t)-(X2(t)-2.0./5.0).^2-(X1(t)-3.0./1.0e1).^2+1.7e1./1.0e1).^(3.0./2.0).*(1.0./2.0)+1.2e1]

depending on time, which should be able to solve like you wanted to do.

Running it as a live script gives:

That's what you wanted to achieve?

Best regards

Stephan

5 个评论
显示 3更早的评论隐藏 3更早的评论

Stephan 2018-5-15

编辑：Stephan 2018-5-15

在 MATLAB Online 中打开

Hi,

the solutions are complex - there seems to be a phase shift - to plot that in you still have to calculate the amplitude and the phase shift - there should be functions for this.

To solve the system, it had to be rewritten as shown:

tspan = 0:0.01:5;
y0 = [0; 0];
[t,X] = ode45(@odefun, tspan, y0)
plot(t, X(1))
hold on
plot(t, X(2))
function ode = odefun(t, X)
ode = zeros(2,1);
ode(1) = (40*X(1) - 12)/(17/10 - X(2) - (X(2) - 2/5)^2 - (X(1) - 3/10)^2 - X(1))^(1/2) - 20*X(1) + ((2*X(1) + 2/5)*(10*X(1)*(2*X(1) - 6/5) + 10*X(2)*(2*X(2) - 6/5)))/(2*(17/10 - X(2) - (X(2) - 2/5)^2 - (X(1) - 3/10)^2 - X(1))^(3/2)) + 12;
ode(2) = (40*X(2) - 12)/(17/10 - X(2) - (X(2) - 2/5)^2 - (X(1) - 3/10)^2 - X(1))^(1/2) - 20*X(2) + ((2*X(2) + 1/5)*(10*X(1)*(2*X(1) - 6/5) + 10*X(2)*(2*X(2) - 6/5)))/(2*(17/10 - X(2) - (X(2) - 2/5)^2 - (X(2) - 3/10)^2 - X(1))^(3/2)) + 12;
end

The solutions are complex - there seems to be a phase shift - to plot that in the plot you still have to calculate the amplitude and the phase shift:

X =
0000 + 0.0000i   0.0000 + 0.0000i
0212 + 0.0000i   0.0214 + 0.0000i
0442 + 0.0000i   0.0448 + 0.0000i
0689 + 0.0000i   0.0704 + 0.0000i
0955 + 0.0000i   0.0983 + 0.0000i
1240 + 0.0000i   0.1287 + 0.0000i
1544 + 0.0000i   0.1616 + 0.0000i
1871 + 0.0000i   0.1975 + 0.0000i
2222 + 0.0000i   0.2367 + 0.0000i
2602 + 0.0000i   0.2798 + 0.0000i
3019 + 0.0000i   0.3280 + 0.0000i
3497 + 0.0000i   0.3842 + 0.0000i
4041 + 0.0000i   0.4493 + 0.0000i
4758 + 0.0000i   0.5358 + 0.0000i
6533 + 0.0339i   0.7352 + 0.0496i
6354 + 0.1948i   0.7471 + 0.2114i
6144 + 0.2632i   0.7507 + 0.3097i
5945 + 0.3075i   0.7453 + 0.3901i
5769 + 0.3378i   0.7337 + 0.4585i
5619 + 0.3584i   0.7175 + 0.5176i
5499 + 0.3721i   0.6983 + 0.5690i
5410 + 0.3807i   0.6770 + 0.6138i
5349 + 0.3855i   0.6540 + 0.6530i
5316 + 0.3874i   0.6300 + 0.6873i
5310 + 0.3872i   0.6052 + 0.7173i
5329 + 0.3855i   0.5802 + 0.7437i
5371 + 0.3827i   0.5548 + 0.7669i
     ...
     ...

I have to sleep now ...

If that helped please accept my anwer in order to help other users finding useful answers on similar problems

Best regards

Stephan

Stephan 2018-5-15

Please note that i have assumed initial conditions:

X1(0)=0

X2(0)=0

You have to check this...i dont have an idea if this is correct.

Best regards

Stephan

Tejas Adsul 2018-5-16

By the way, how did you write down the expression of ode(1) and ode(2) in the function odefun?

请先登录，再进行评论。

Answer 2

Stephan 2018-5-15

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/400867-help-solving-a-system-of-differential-equations#answer_320338

编辑：Stephan 2018-5-15

在 MATLAB Online 中打开

Hi,

did i understand right:

x1=0.3;
x3=0.9; 
y1=0.4;
y3=0.8;
A=10; 
p=10;
syms X1 X2;
num = ((x3-X1)-(X1-x1))*X1 + ((y3-X2)-(X2-y1))*X2;
den = sqrt((x3-X1)-(X1-x1)^2 + (y3-X2)-(X2-y1)^2);
Em = -A*(X1-x1)*(x3-X1) - A*(X2-y1)*(y3-X2) + p*num/den;
dEm_dX1 = diff(Em,X1)
dEm_dX2 = diff(Em,X2)
dX1dt = -dEm_dX1
dX2dt = -dEm_dX2

This is what you wanted to do?

gives:

dX1dt =
(40*X1 - 12)/(17/10 - X2 - (X2 - 2/5)^2 - (X1 - 3/10)^2 - X1)^(1/2) - 20*X1 + ((2*X1 + 2/5)*(10*X1*(2*X1 - 6/5) + 10*X2*(2*X2 - 6/5)))/(2*(17/10 - X2 - (X2 - 2/5)^2 - (X1 - 3/10)^2 - X1)^(3/2)) + 12
dX2dt =
(40*X2 - 12)/(17/10 - X2 - (X2 - 2/5)^2 - (X1 - 3/10)^2 - X1)^(1/2) - 20*X2 + ((2*X2 + 1/5)*(10*X1*(2*X1 - 6/5) + 10*X2*(2*X2 - 6/5)))/(2*(17/10 - X2 - (X2 - 2/5)^2 - (X1 - 3/10)^2 - X1)^(3/2)) + 12

Can this be correct?

Best regards

Stephan

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Tejas Adsul 2018-5-15

Yes, but this is the first part, which I have been able to do. The second part involves actually solving the ODE numerically (I don't think it can be solved analytically) and getting X1 and X2 as functions of time. In the end, I would like to plot (t,X1) and (t,X2).

请先登录，再进行评论。