6 views (last 30 days)

function df = myfcn(t,f)

df = zeros(3,1);

df(1) = f(2);

df(2) = f(3) ;

df(3) = (-.5*f)*f(3) ;

tspan = [0 10];

f0 = [0 0];

f10= [10 1];

[t,f] = ode45(@myfcn, tspan, f0, f10);

plot(t,f,'-o',t,f(:,1),'-.',t,f(:,2),'-o')

end

Stephan
on 5 Dec 2019

Edited: Stephan
on 5 Dec 2019

There are at least 4 problems in your code:

1. Split your main code and your ode function:

tspan = [0 10];

f0 = [0 0];

f10= [10 1];

[t,f] = ode45(@myfcn, tspan, f0, f10);

plot(t,f,'-o',t,f(:,1),'-.',t,f(:,2),'-o')

function df = myfcn(t,f)

df = zeros(3,1);

df(1) = f(2);

df(2) = f(3) ;

df(3) = (-.5*f).*f(3) ;

end

2. There is an initial condition needed for every ode. Since you have 3 equations you have to give f0 as 1x3 vector, not as 1x2.

3. If you want to use f10 as boundary condition ode45 is the wrong tool. To solve boundary value problems use bvp4c. For the usage of ode45 leave it away.

4. The last equation has f*... - f is a 3x1 vector, which will cause problems, because df(3) is expected to be a scalar. Maybe you mean f(1) for example.

Making these corrections should give you an error free code. If it is meaningful in the end i can not say.

Stephan
on 5 Dec 2019

Stephan
on 6 Dec 2019

Your code is ok - just the initial condition is not useful. With trying a little bit around you can find 1/3 as near to the truth, if we want to meet the boundary condition at t=10:

f0 = [0 0 1/3];

But i think your instructor makes something wrong when he tells you to use ode45 for this. As already stated this is a boundary value problem, so the right tool is bvp4c to tackle this. If there were only initial values ode45 would be the right choice.

However, once we tried out that 1/3 is a good estimation for the initial conditions when using ode45, we can compare both approaches:

% Using ode 45 needs a very good guess of the initial

% value for d2f - this one is near to good but not perfect

tspan = [0 10];

f0 = [0 0 1/3];

[t,f] = ode45(@myfcn, tspan, f0);

result_of_boundary_condition_ode45 = f(end,2)

% bvp4c allows to set the boundary conditions as given,

% there is no need to try out for a good initial condition

% like in ode45

xmesh = linspace(0, 10, 11);

solinit = bvpinit(xmesh, [0 0 0]);

sol = bvp4c(@myfcn, @bcfcn, solinit);

result_of_boundary_condition_bvp4c = sol.y(2,end)

% plot results

plot(t,f(:,1),'b',t,f(:,2),'r',t,f(:,3),'g')

hold on

plot(sol.x,sol.y(1,:),'ob',...

sol.x,sol.y(2,:),'or',...

sol.x,sol.y(3,:),'og')

hold off

legend('f - ode45','df - ode45','d2f - ode45',...

'f - bvp4c','df - bvp4c','d2f - bvp4c',...

'location','Northwest')

title('Results of ode45 and bvp4c')

xlabel('$\eta$','Interpreter','latex','FontSize',15);

ylabel('F, dF, d2F')

% ode function (same for both approaches)

function df = myfcn(~,f)

df = zeros(3,1);

df(1) = f(2);

df(2) = f(3) ;

df(3) = (-.5*f(1))*f(3) ;

end

% Boundary conditions (needed for bvp4c)

function res = bcfcn(ya,yb)

res = [

ya(1)

ya(2)

yb(2)-1];

end

This looks quiet good, but it is not perfect. Note what happens if you change the initial condition to lets say 0.5:

f0 = [0 0 0.5];

If you look at the plot now, you know why bvp4c should be used.

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.