Second Order ODE solved with ODE45

Hi!

I am struggling with the task to solve a 2'nd order Ode with ODE45. Please help as I'm struggling with this.

Equation as given in the task:

y'' + pi*y^(x/3)*(2y' sin pi*x + pi*y cos pi*x) − y/9 = 0, y(0) = and y'(0) = -⅓

My code:

%Lab 2 uppgift 3
clear all
close all
clc

[t,y] = ode45(@odefun, t, [0, 2.5], [1, -1/3]);

plot(t,y(:,1),'-',t,y(:,2),'--')

title('Title');

xlabel('time t');

ylabel('solution y');

legend('y_1','y_2')

%The function called odefun.

function [dydx] = odefun(t,y)

%x = (0:0.1:2.5);

dydt=[y(2); pi*y(1).*exp(x/3)*2*y(2)*sin(pi*x) + x*y(1)*cos(pi*x) - (y(1)./9)x];

end

Thanks! /Henrik

Answer 1

[Star Strider]

Your equation implies (to me) that y is a function of x, not t. So your function should probably be:

function [dydx] = odefun(x,y)

Also, there is an error in this line at the very end:

dydt=[y(2); pi*y(1).*exp(x/3)*2*y(2)*sin(pi*x) + x*y(1)*cos(pi*x) - (y(1)./9)x];

You might want:

dydt=[y(2); pi*y(1).*exp(x/3)*2*y(2)*sin(pi*x) + x*y(1)*cos(pi*x) - (y(1)./9).*x];

or some other operator, since MATLAB will throw an error without an operator there.

When I ran it (with the change in the function line, vectorising, and putting a multiplication operator before the x at the end, it ran.

[Henrik]

Thanks a lot star strider!

The plot I get doesn't look right but I get it running at least. However with an error message saying:

"Warning: Failure at t=1.425454e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552714e-15) at time t."

If you have any ideas to get rid of it I'm even more grateful.

Thanks //Henrik

[Star Strider]

My pleasure!

I re-coded your equation and found some errors. This seems to produce a much more reasonable result:

odefun = @(x,y) [y(2); -pi*y(1).^(x/3).*2.*y(2).*sin(pi*x) - x.*y(1).*cos(pi*x) + (y(1)./9).*x];

[t,y] = ode15s(odefun, [0, 2.5], [1, -1/3]);

plot(t,y(:,1),'-',t,y(:,2),'--')  
title('Title');  
xlabel('time t');  
ylabel('solution y');  
legend('y_1','y_2','Location','NW')

I also repositioned your legend so that it would not be on top of your plotted curves.