solve and plot a system of nonlinear 2nd order differential equations

Question

Zhen Zhen 2019-3-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/452333-solve-and-plot-a-system-of-nonlinear-2nd-order-differential-equations

编辑： Lewis Fer 2021-6-12

采纳的回答： Teja Muppirala

在 MATLAB Online 中打开

hi there,

I'm trying to plot a graph of

against

with the following equations of motion:

I've tried dsolve and ode45 yet there always seems to be some problems. I think ode45 might work better because apparently it would be easier to plot the graph by using some numerical method?

Here's my failed attempt to solve it: (I just set some variables to be 1 to make the problem easier here)

syms theta(t) phi(t) psi(t) C
dtheta = diff(theta , t);
d2theta = diff(theta , t , 2);
dphi = diff(phi , t);
%dpsi = diff(psi , t);
%alpha = C * (dpsi + dphi * cos(theta));
%beta = alpha * cos(theta) + dphi * (sin(theta))^2;
alpha = 1;
beta = 1;
dpsi = 1;
eqn1 = dphi == (beta - alpha * cos(theta)) / (sin(theta))^2 ; %equations of motion
eqn2 = d2theta == (dphi*(dphi * cos(theta) - alpha)+1)*sin(theta) ; %equations of motion
eqns = [eqn1 , eqn2];
%cond = [dpsi == 1];
[thetaSol(t) phiSol(t)] = dsolve(eqns)

Thanks a lot for your help and time in advance!

Cheers,

Jane

2 个评论
显示无隐藏无

Walter Roberson 2019-3-26

在 MATLAB Online 中打开

There are at least four solutions to that.

Two of them are pretty abstract, involving integrals of roots of an equation -- not closed form solutions but rather a description of what properties the solution function would have to have.

The other two solutions are closed form:

theta(t) = arctan(sqrt(2*diff(phi(t), t) - 1)/diff(phi(t), t), (-diff(phi(t), t) + 1)/diff(phi(t), t))

and the same except the negative of the first parameter.

MATLAB is not powerful enough to arrive at these solutions.

Walter Roberson 2019-3-26

Have a look at odeFunction() and in particular the first example. It shows you the functions you need to use in order to convert the symbolic forms into something you can call with ode45.

Use the options structure created by odeset to designate an OutputFcn . You might want to use @odephas2 to construct a 2D phase plot, perhaps.

请先登录，再进行评论。

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

Answer 1

Teja Muppirala 2019-3-26

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/452333-solve-and-plot-a-system-of-nonlinear-2nd-order-differential-equations#answer_367408

编辑：Teja Muppirala 2019-3-26

在 MATLAB Online 中打开

If you just need a plot and not a closed-form solution, then I'd recommend just using ODE45 without worrying about symbolic stuff. This is an example of how to solve this using ODE45 for initial conditions psi(0) = 0, theta(0) = 0, thetadot(0) = 1 over the time span [0 10].

function doODE
a = 1; % alpha
b = 1; % beta
% Define variables as:
% Y(1): phi
% Y(2): theta
% Y(3): thetadot
ic = [0;0;1]; % Pick some initial Conditions 
tspan = [0 10];
[tout, yout] = ode45(@deriv,tspan,ic);
subplot(1,3,1), plot(tout,yout(:,1));  title('psi(t)')
subplot(1,3,2), plot(tout,yout(:,2)); title('theta(t)')
subplot(1,3,3), plot(yout(:,1),yout(:,2)); title('theta vs psi')
    function dY = deriv(t,Y)
        
        dY = [dpsi(Y(2)); ...
            Y(3); ...
            [dpsi(Y(2))*(dpsi(Y(2))*cos(Y(2))-a)+1]*sin(Y(2)) ];
    end
    function dp = dpsi(th)
        
        dp = (b - a*cos(th))./(sin(th)).^2 ;
        if isnan(dp) % Need to take care at theta = 0
            dp = 0.5;
        end
    end
end

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

Zhen Zhen 2019-3-28

Thank you very much! it worked well!

I was wondering, if now I need to use the real alpha and beta (alpha = C * (dpsi + dphi * cos(theta)); beta = alpha * cos(theta) + dphi * (sin(theta))^2;) I know they are constants. Where should I put them in the code?

And if I need to print one more constant (say, E = 0.5 * (Y(3))^2 + (Y(1))^2 * (sin(Y(2)))^2) +0.5* a^2 + cos(Y(2)) ) where should I put it?

Thanks so much for your help!

请先登录，再进行评论。