Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
part c
  • Loading branch information
nin13001 committed Dec 6, 2017
1 parent fdd7e99 commit 0908032
Show file tree
Hide file tree
Showing 5 changed files with 89 additions and 0 deletions.
22 changes: 22 additions & 0 deletions 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
22 changes: 22 additions & 0 deletions 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
28 changes: 28 additions & 0 deletions README.md
Expand Up @@ -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)
Binary file added assets/README-23f54e97.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
17 changes: 17 additions & 0 deletions 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')

0 comments on commit 0908032

Please sign in to comment.