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

1 次查看（过去 30 天）
Kaylene Widdoes
Kaylene Widdoes 2016-3-11
评论： Star Strider 2016-3-11
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 个评论
Kaylene Widdoes
Kaylene Widdoes 2016-3-11
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, ... >>

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2016-3-11
编辑：Star Strider 2016-3-11
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 个评论
Kaylene Widdoes
Kaylene Widdoes 2016-3-11
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 2016-3-11
Thank you!
Yes. I switched it to ode15s and it gave a much better solution. (Sorry for the delay — life intrudes.)

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by