diff --git a/Prob2ptC.asv b/Prob2ptC.asv new file mode 100644 index 0000000..83b6c4a --- /dev/null +++ b/Prob2ptC.asv @@ -0,0 +1,22 @@ +function y=Prob2ptC(ode,dt,y0,tlimits) +maxit=100; +t=[tlimits(1):dt:tlimits(2)]; +y=zeros(length(t),1); +y(1)=y0; + for i=2:length(t) + dy=ode(t(i-1),y(i-1)); + y(i)=y(i-1)+dt*dy; + dyc=(dy+ode(t(i),y(i)))/2; + y(i)=y(i-1)+dt*dyc; + nits=1; + + while(1) + nits=nits+1; + yold=y(i); + dyc=(dy+ode(t(i),y(i)))/2; + y(i)=y(i-1)+dt*dyc; + ea=abs(yold-y(i))/y(i)*100; + if ea<.0001|nits>maxit;break;end; + end; + end +end \ No newline at end of file diff --git a/Prob2ptC.m b/Prob2ptC.m new file mode 100644 index 0000000..83b6c4a --- /dev/null +++ b/Prob2ptC.m @@ -0,0 +1,22 @@ +function y=Prob2ptC(ode,dt,y0,tlimits) +maxit=100; +t=[tlimits(1):dt:tlimits(2)]; +y=zeros(length(t),1); +y(1)=y0; + for i=2:length(t) + dy=ode(t(i-1),y(i-1)); + y(i)=y(i-1)+dt*dy; + dyc=(dy+ode(t(i),y(i)))/2; + y(i)=y(i-1)+dt*dyc; + nits=1; + + while(1) + nits=nits+1; + yold=y(i); + dyc=(dy+ode(t(i),y(i)))/2; + y(i)=y(i-1)+dt*dyc; + ea=abs(yold-y(i))/y(i)*100; + if ea<.0001|nits>maxit;break;end; + end; + end +end \ No newline at end of file diff --git a/README.md b/README.md index 3ce5658..a5c0450 100644 --- a/README.md +++ b/README.md @@ -19,3 +19,31 @@ end plot(x,exp(-x),'k-',x,y_n,'o'); ``` ![](assets/README-66741291.png) + +## Part C + +``` +function y=Prob2ptC(ode,dt,y0,tlimits) +maxit=100; +t=[tlimits(1):dt:tlimits(2)]; +y=zeros(length(t),1); +y(1)=y0; + for i=2:length(t) + dy=ode(t(i-1),y(i-1)); + y(i)=y(i-1)+dt*dy; + dyc=(dy+ode(t(i),y(i)))/2; + y(i)=y(i-1)+dt*dyc; + nits=1; + + while(1) + nits=nits+1; + yold=y(i); + dyc=(dy+ode(t(i),y(i)))/2; + y(i)=y(i-1)+dt*dyc; + ea=abs(yold-y(i))/y(i)*100; + if ea<.0001|nits>maxit;break;end; + end; + end +end +``` +![](assets/README-23f54e97.png) diff --git a/assets/README-23f54e97.png b/assets/README-23f54e97.png new file mode 100644 index 0000000..a16cda4 Binary files /dev/null and b/assets/README-23f54e97.png differ diff --git a/prob2plot.m b/prob2plot.m new file mode 100644 index 0000000..be9521c --- /dev/null +++ b/prob2plot.m @@ -0,0 +1,17 @@ +set(0, 'defaultAxesFontSize', 16) +set(0,'defaultTextFontSize',14) +set(0,'defaultLineLineWidth',3) +dt = 0.1; +t = [0:dt:3]; +y = zeros(size(t)); +y(1) = 1; +for i = 2:length(t) + dy = prob2_ode(t,y(i-1)); + y(i) = y(i-1)+dt*dy; +end +a = figure(1); +y_heun = Prob2ptC (@(t,y) prob2_ode(t,y), 0.1, 1, [0 3]); +plot(t,exp(-t),'k-',t,y,'Xg',t,y_heun,'Or'); +legend('Analytical','Euler','Heun','Location','Northeast') +xlabel('Time') +ylabel('Position') \ No newline at end of file