2nd Order Non-Linear Equation Numerical Solution and Plot

Hi All,

I am looking to numberically solve and graph the following equation:

With initial conditions:

x=x(t); A and B are constant coefficients. I've been trying to use ODE45 but can't seem to get it to work properly.

Any help is appreciated!

[William Rose]

@Conlen O'Brien, can you share what you have tried that didn;t work?

[Conlen O'Brien]

@William Rose

I may have just figured it out, see below:

function dxdt= odefcn(t,x,A,B)
    dxdt= zeros(2,1);
    dxdt(1)=x(2);
    dxdt(2)=-A+B*x(2)^2
end
--------------------------------------
clear
clc

A=9.81;
B=3;
tspan = [0 10];
x0=[-25 0];
[t,x]=ode45(@(t,x) odefcn(t,x,A,B),tspan,x0);
plot(t,x(:,1),'-o',t,x(:,2),'-.')

Note that the function is saved in a seperate file. It gives me the following plot:

Am I correct thinking that the blue plot is x(t), and the orange plot is x'(t)?

[William Rose]

@Conlen O'Brien, Good job!

[William Rose]

@Conlen O'Brien,

Here are a couple of commands to plot results:

A=9.81; B=3;
tspan=[0 10];
x0=[-25 0];
[t,x] = ode45(@(t,x) [x(2);-A+B*x(2).^2], tspan, x0);

Plot results

subplot(211); plot(t,x(:,1),'-r.'); ylabel('x(t)'); grid on
subplot(212); plot(t,x(:,2),'-r.'); ylabel('dx/dt'); grid on; xlabel('Time')

You can see, by analyzing the original differential equation, that a "steady state" solution is

dx/dt=sqrt(A/B). Since dx/dt=constant, d2x/dt2=0, so this satisfies the original differential equation. Does the numerical solution above agree with this? Yes. sqrt(A/B)=sqrt(9.81/3)=+1.81 or -1.81. The plotted solution agrees with this: dx/dt=-1.81 in the "steady state".

Answer 1

[William Rose]

@Conlen O'Brien,

tspan=[0 10];  x0=[-25 0];
A=1; B=1;
[t,x] = ode45(@(t,x) [x(2);-A+B*x(2).^2], tspan, x0);

Do you think the results make sense?

[Sam Chak]

Hi @Conlen O'Brien

You can assign a display name to each state and toggle the legends. Also, if A and B are fixed constants, then place them inside the odefcn() function so that you don't need to assign the values to A and B every time you call the ode45 solver.

tspan = [0 10];
x0    = [-25 0];
[t,x] = ode45(@odefcn, tspan, x0);

plot(t, x(:,1), 'blue', 'DisplayName', 'x(t)'),  hold on
plot(t, x(:,2), 'red', 'DisplayName', 'x''(t)'), grid on
xlabel('t'), ylabel('\bf{x}')
legend show

function dxdt= odefcn(t, x)
    A       = 9.81;
    B       = 3;
    dxdt    = zeros(2,1);
    dxdt(1) = x(2);
    dxdt(2) = - A + B*x(2)^2;
end