Problem 2 and 3

HW6Problem2.m

@@ -0,0 +1,20 @@
+%Part b.
+dt = [0.001];
+figure()
+hold on
+for k = 1:length(dt)
+    t = [0:dt:3]';
+    y_n = zeros(size(t));
+    y_n(1) = 1; %initial condition
+    for i = 2:length(t)
+        dy = ode2(t, y_n(i-1));
+        y_n(i) = y_n(i-1) + dt(k)*dy;
+    end
+
+    plot(t, exp(-t), 'k-', t, y_n, 'o-');
+end
+
+%Part c.
+figure()
+y_heun = heun_sol_order1(@(t,y) ode2(t,y), 1, 1, [0 3]);
+plot([0:1:3], y_heun);

HW6Problem3.m

@@ -0,0 +1,23 @@
+%Part b.
+dy = ode3(0, [1;0]);
+dt = [0.1 0.01 0.001];
+figure()
+hold on
+for k = 1:length(dt)
+    t = [0:dt(k):3]';
+    y_n = zeros(length(t), 2);
+    y_n(1, :) = [1 0]; %initial condition
+    for i = 2:length(t)
+        dy = ode3(t, y_n(i-1, :));
+        y_n(i, :) = y_n(i-1, :) + dt(k)*dy;
+    end
+
+    plot(t, cos(3*t), 'k-', t, y_n(:,1), 'o-');
+end
+
+%Part c.
+figure()
+dt = 0.001;
+t = [0:dt:3];
+y_heun = heun_sol_order2(@(t,y) ode3(t,y), dt, [1 0], [0 3]);
+plot(t, cos(3*t), 'k-', t, y_heun(:,1), 'o-');

heun_sol.m

@@ -0,0 +1,27 @@
+function y = heun_sol_order1(ode, dt, y0, tlimits)
+    %Solves ODE using heun's method
+    %returns y = solution to ODE
+    %ode = differential equation to be solved
+    %dt = size of time step
+    %y0 = initial conditions
+    %tlimits = upper and lower bounds of desired time interval
+    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;
+        n = 0;
+        error = 100;
+        while(1)
+            n = n+1;
+            yOld = y(i);
+            dyc = (dy + ode(t(i), y(i))) / 2;
+            y(i) = y(i-1) + dt*dyc;
+            error = abs(y(i)-yOld)/y(i)*100;
+            if error <= 0.00001 | n>100, break, end
+        end
+    end
+end

heun_sol.m~

@@ -0,0 +1,23 @@
+function y = heun_sol(ode, dt, y0, tlimits)
+    %Solves ODE using heun's method
+    %
+    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;
+        n = 0;
+        error = 100;
+        while(1)
+            n = n+1;
+            yOld = y(i);
+            dyc = (dy + ode(t(i), y(i))) / 2;
+            y(i) = y(i-1) + dt*dyc;
+            error = abs(y(i)-yOld)/y(i)*100;
+            if error <= 0.00001 | n>100, break, end
+        end
+    end
+end

heun_sol_order1.m

@@ -0,0 +1,27 @@
+function y = heun_sol_order1(ode, dt, y0, tlimits)
+    %Solves ODE using heun's method
+    %returns y = solution to ODE
+    %ode = differential equation to be solved
+    %dt = size of time step
+    %y0 = initial conditions
+    %tlimits = upper and lower bounds of desired time interval
+    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;
+        n = 0;
+        error = 100;
+        while(1)
+            n = n+1;
+            yOld = y(i);
+            dyc = (dy + ode(t(i), y(i))) / 2;
+            y(i) = y(i-1) + dt*dyc;
+            error = abs(y(i)-yOld)/y(i)*100;
+            if error <= 0.00001 | n>100, break, end
+        end
+    end
+end

heun_sol_order2.m

@@ -0,0 +1,27 @@
+function y = heun_sol_order2(ode, dt, y0, tlimits)
+    %Solves ODE using heun's method
+    %returns y = solution to ODE
+    %ode = differential equation to be solved
+    %dt = size of time step
+    %y0 = initial conditions
+    %tlimits = upper and lower bounds of desired time interval
+    t = [tlimits(1):dt:tlimits(2)];
+    y = zeros(length(t), length(y0));
+    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;
+        n = 0;
+        error = 100;
+        while(1)
+            n = n+1;
+            yOld = y(i, :);
+            dyc = (dy + ode(t(i), y(i, :))) / 2;
+            y(i, :) = y(i-1, :) + dt*dyc;
+            error = abs(y(i, :)-yOld)./y(i, :)*100;
+            if max(error) <= 0.00001 | n>100, break, end
+        end
+    end
+end

ode2.m

@@ -0,0 +1,4 @@
+function dy = ode2(t,y)
+    %differential equation from hw6 problem 2
+    dy = -y;
+end

ode3.m

@@ -0,0 +1,6 @@
+function dy = ode3(t,y)
+    %differential equation from hw6 problem 3
+    dy = zeros(size(y));
+    dy(1) = y(2);
+    dy(2) = -9*y(1);
+end