diff --git a/HW6Problem2.m b/HW6Problem2.m new file mode 100644 index 0000000..c186f79 --- /dev/null +++ b/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); diff --git a/HW6Problem3.m b/HW6Problem3.m new file mode 100644 index 0000000..ecd77ce --- /dev/null +++ b/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-'); \ No newline at end of file diff --git a/heun_sol.m b/heun_sol.m new file mode 100644 index 0000000..daa8961 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/heun_sol.m~ b/heun_sol.m~ new file mode 100644 index 0000000..71a2d2e --- /dev/null +++ b/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 \ No newline at end of file diff --git a/heun_sol_order1.m b/heun_sol_order1.m new file mode 100644 index 0000000..daa8961 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/heun_sol_order2.m b/heun_sol_order2.m new file mode 100644 index 0000000..25f9f63 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/ode2.m b/ode2.m new file mode 100644 index 0000000..98e0989 --- /dev/null +++ b/ode2.m @@ -0,0 +1,4 @@ +function dy = ode2(t,y) + %differential equation from hw6 problem 2 + dy = -y; +end \ No newline at end of file diff --git a/ode3.m b/ode3.m new file mode 100644 index 0000000..832cafa --- /dev/null +++ b/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 \ No newline at end of file