ODE45 - How do I set up the conditions so that I'm not getting an error on 0?

1 vue (au cours des 30 derniers jours)
I have to choose an appropriate ode solver for y'(t) = te^(-1000*t) - 10*y(t) with y(0) = 0 for 0<=t<=0.5. I also have to choose appropriate values for the absolute and relative tolerance to get a solution. I am getting an error with the code I have tried to start running to plot a vector solution of it. Also, can anyone confirm that the ODE45 solver and tolerance values are the correct choice? Thanks!
a = 0;
b = .5;
y(1) = 0;
f1 = @(t,y) (t*exp(-1000*t)) - (10*y(t));
options = odeset('RelTol',10^-7,'AbsTol',10^-9);
z1 = cputime;
[T Y] = ode45(f1,[a b],0,options);
z2 = cputime;
fprintf('ODE 45 CPU time: %f seconds. Number of steps: %d\n',z2-z1,length(T))
plot(T,Y,'r')
  2 commentaires
Steven Lord
Steven Lord le 11 Mar 2016
What is the full text of the error you received?
Kaylene Widdoes
Kaylene Widdoes le 11 Mar 2016
Attempted to access y(0); index must be a positive integer or logical. Error in @(t,y)(t*exp(-1000*t))-(10*y(t))
Error in odearguments (line 87) f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0. Error in ode45 (line 113) [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ... >>

Connectez-vous pour commenter.

Réponse acceptée

Star Strider
Star Strider le 11 Mar 2016
Modifié(e) : Star Strider le 11 Mar 2016
This runs for me without error:
a = 0;
b = .5;
y(1) = 0;
f1 = @(t,y) (t.*exp(-1000*t)) - (10*y);
options = odeset('RelTol',10^-7,'AbsTol',10^-9);
z1 = cputime;
[T Y] = ode45(f1,[a b],0,options);
z2 = cputime;
fprintf('ODE 45 CPU time: %f seconds. Number of steps: %d\n',z2-z1,length(T))
plot(T,Y,'r')
The only changes I made were to replace ‘y(t)’ with ‘y’ in ‘f1’, and vectorise it.
  2 commentaires
Kaylene Widdoes
Kaylene Widdoes le 11 Mar 2016
Awesome! Ok - so that looks like a really rigid graph - would I be better off switching to a different solver or tolerance values to accommodate for the stiffness?
Star Strider
Star Strider le 11 Mar 2016
Thank you!
Yes. I switched it to ode15s and it gave a much better solution. (Sorry for the delay — life intrudes.)

Connectez-vous pour commenter.

Plus de réponses (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by