Blank figures appearing for plots

2 views (last 30 days)
Phoebe Tyson
Phoebe Tyson on 10 Mar 2020
Commented: Image Analyst on 10 Mar 2020
I've create an adaptive time-stepping function and I'm trying to plot the graphs but it is just creating blank figures.
Here is my function:
function [tvec,yvec,dtvec]=adaptive_euler(f,ytrue,y0,t0,dt0,dtmin,dtmax,tmax,alpha,beta,gamma,TOL)
y=zeros(10^5,1);
y(1)=y0;
dt=zeros(10^5,1);
t=zeros(10^5,1);
dt(1)=dt0;
yt=y0;
yt1=y0;
t(1)=t0;
n=1;
while (t(n)<tmax && dt(n)<dtmax && dt(n)>dtmin)
yt=yt+dt(n)*f(t(n),yt);
yt1=yt1+(dt(n)/2)*f(t(n),yt1);
yt2=yt1+(dt(n)/2)*f(t(n),yt1);
Et2=(yt2-yt)/(dt(n));
c=TOL/norm(Et2);
if norm(Et2) < TOL
y(n)=yt2;
t(n+1)=t(n)+dt(n);
n=n+1;
end
dt(n) = min(max(alpha*c*dt(n-1),gamma*dt(n-1)),beta*dt(n-1));
end
tvec=zeros(n,1);
yvec=zeros(n,1);
dtvec=zeros(n,1);
for i=1:n
tvec(i)=t(i);
yvec(i)=y(i);
dtvec(i)=dt(i);
end
yreal=zeros(n,1);
error=zeros(n,1);
for i=1:n
yreal(i)=ytrue(tvec(i));
error(i)=abs(yvec(i)-yreal(i));
end
T1=table(n,error(n));
figure(1)
plot(tvec,yvec,'.',tvec,yreal,'--');
legend('Approx:','True:');
title('Graph of the Adaptive Eulers Method approximation of dy/dx=-exp(2*t)*y^2 against real values');
figure(2)
plot(tvec,dtvec)
title('Graph of the change in dt over t');
end
And here is the script:
f=@(t,y)-exp(2*t)*y^2;
t0=0;
y0=2;
dt0=0001;
dtmin=10^-10;
dtmax=0.5;
tmax=10;
alpha=0.84;
beta=4;
gamma=0.1;
TOL=0.1;
ytrue=@(t)2*exp(-2*t);
[tvec,yvec,dtvec]=adaptive_euler(f,ytrue,y0,t0,dt0,dtmin,dtmax,tmax,alpha,beta,gamma,TOL)
Thank you!
  4 Comments
Adam
Adam on 10 Mar 2020
But you set them up as vectors of length n and then you copy from 1 to n from another vector into these.
Phoebe Tyson
Phoebe Tyson on 10 Mar 2020
Because in the original vectors, they are much longer and everything after n is a zero, n is the number of iterations my while loop runs for, while the length of the original vectors is much larger

Sign in to comment.

Answers (1)

Image Analyst
Image Analyst on 10 Mar 2020
You never enter the while loop
while (t(n)<tmax && dt(n)<dtmax && dt(n)>dtmin)
because dt(1) is 1 and dtmax is 0.5. Pass in a different dtmax.
  2 Comments
Phoebe Tyson
Phoebe Tyson on 10 Mar 2020
Ah, it was supposed to be 0.001 instead of 1. now there seems to be an issue with
dt(n) = min(max(alpha*c*dt(n-1),gamma*dt(n-1)),beta*dt(n-1));
Saying "Array indices must be positive integers or logical values."
Image Analyst
Image Analyst on 10 Mar 2020
1 is also not less than 0.001, so that won't enter the loop either. So what did you pass in?

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!