MATLAB Answers

I am trying to use ODE45, but I keep getting an error with myfcn saying I don't have enough inputs?

2 views (last 30 days)
Hannah Middlestead
Hannah Middlestead on 5 Dec 2019
Edited: Stephan on 6 Dec 2019
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

  0 Comments

Sign in to comment.

Answers (1)

Stephan
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.

  4 Comments

Show 1 older comment
Stephan
Stephan on 5 Dec 2019
You dont have to give up. I think we can solve your problem. Post what you have tried, which errors you get and it would be good to see the mathematical model in latex form to get an idea of your system.
Hannah Middlestead
Hannah Middlestead on 5 Dec 2019
Alrighty...
This is a picture of problem (don't have Latex):HW5.JPG
This is what I tried after your suggestions:
tspan = [0 10];
f0 = [0 0 0];
f1=[1 0 0];
f10= [10 0 1];
[t,f] = ode45(@myfcn, tspan, f0,f1);
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(1))*f(3) ;
end
Stephan
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.

Tags