MATLAB: How to solve this
differential equationseuler's methodhomework
Related Solutions
tspan = [0 4];y0 = 2;[t1,y1] = ode45(@(t,y) 4*exp(0.8*t)-0.5*y,tspan,y0,h);figure('name','Reda','color','w');plot(t1,y1, 'r-','LineWidth',2);grid on;xlabel('Time');ylabel('Solution y(Time)');title('And you instead of a blanket substitute solver Euler');
The first thing you need to do is to write the ode as two first order equations and use the code below. You will be required to supply two initial conditions for the 1s order equations. Use the one that you are given plus another of your choice.
function [t,x,y,N] = Runge2_2eqs(f1,f2,to,tfinal,xo,yo,h) % This function implements the Rk2 method.
t = to; N = ceil((tfinal-to)/h); x = zeros(1,N); y = zeros(1,N) ; x(1) = xo; y(1) = yo; for i = 1:N t(i+1) = t(i)+h; Sx1 = f1(t(i),x(i),y(i)); Sy1 = f2(t(i),x(i),y(i)); Sx2 = f1(t(i)+h, x(i)+Sx1*h, y(i)+Sy1*h); Sy2 = f2(t(i)+h, x(i)+Sx1*h, y(i)+Sy1*h); x(i+1) = x(i) + h/2*(Sx1+Sx2); y(i+1) = y(i) + h/2*(Sy1+Sy2); endend
This is the mfile.
xo = 1;yo = 0;h = [.1 0.03];to = 0;tfinal = 20;M = ceil((tfinal-to)/h(2));dx1 = @(t,x1,x2) x2;dx2 = @(t,x1,x2) 6*x2 -13*x1 + t + 3*sin(t);% When you reduce the equation to two first order, x will be the solution
% of the ode, i.e x'' and y is its derivative, x'.
for i = 1: length(h)if (i== 1) % This for the case when h = 0.1
[t,x,y,N] = Runge2_2eqs(dx1,dx2,to,tfinal,xo,yo,h(i)); y1 = x; y2 = y;else % and for the case when h = 0.03
[t,x,y,N] = Runge2_2eqs(dx1,dx2,to,tfinal,xo,yo,h(1)); x3 = x; x4 = y;endendt1 = t(1):(t(end)-t(1))/(M-1):t(end);figure(1);subplot(121)plot(t1,y1, '-o')hold onplot(t1,y2,'-o')legend('Dfx1','Dfx2')title('Solution to two systems of ODEs using RK2, h= 0.1')xlabel('x')ylabel('y')xlim([to tfinal])gridsubplot(122)plot(t,x3,'linewidth',2,'color','b')hold onplot(t,x4,'linewidth',2,'color','r')legend('Dfx1','Dfx2')title('Solution to two systems of ODEs using RK2, h = 0.03')xlabel('x')ylabel('y')xlim([to tfinal])grid% Using ode 45 just to prove that the solution with RK2 is correct.
F = @(t,y) [ y(2); (6*y(2) -13*y(1) + t + 3*sin(t)) ];t0 = 0;tf = 20;delta = (tf-t0)/(201-1); tspan = t0:delta:tf;ic = [1 0];[t,y] = ode45(F, tspan, ic);figureplot(t,y(:,1),'-o')hold onplot(t,y(:,2),'-o')a = title('Using ode45');legend('x','x_{prime}');set(a,'fontsize',14);a = ylabel('y');set(a,'Fontsize',14);a = xlabel('t [0 20]');set(a,'Fontsize',14);xlim([t0 tf])gridgrid minor;
Best Answer