# Sine equation in Euler Method ODE

2 views (last 30 days)
Chris Horne on 31 Mar 2022
Commented: Chris Horne on 1 Apr 2022
I am trying to code the following equation
y = SQRT(1+x)*cos(x^2)
The exact solution is so much different from the Euler result. I am njot sure if it is axes or code typo. Any help is appreciated. Thanks
% 2D SYSTEM SIN-COS EQUATIONS USING FORWARD EULER METHOD
% % Initial conditions and setup
% x(t)= sqrt(t+1)*cos(t).^2
t(1)=0; % initial t
x(1)=1; y(1)=0; ; % initial x,y
dt=0.005; % time step
nn=10000; % number of time steps in radians
for i=1:nn-1 % loop counter
dx=(sqrt(i+1)*cos(i.^2))./2*(i+1) - 2*i*(sqrt(i+1)*sin(i.^2)); %x' equation
dy=(sqrt(i+1)*sin(i.^2))./2*(i+1) - 2*i*(sqrt(i+1)*cos(i.^2)); %y' equation
x(i+1)=x(i)+dt*dx; % Find new x
y(i+1)=y(i)+dt*dy; % Find new y
t(i+1)=t(i)+dt; % Find new t
end
% exact solution without Euler
tt = linspace(0,10000,length(x));%set the t,x vectors to same lengthtt=0:0.1:10000;
x_exact = sqrt(1+tt).*cos(tt.^2); %????? typo?
% Plot x vs. t with exact solution
figure
t = linspace(0,10000,length(x)); %set the t,x vectors to same length
plot(t,x,'b')
%axis([0 10000 -5000 5000])
hold on
title('2D System sqrt(1+t)cos equation vs exact soln')
xlabel('t')
ylabel('x')
plot(t,x_exact,'r')
##### 2 CommentsShowHide 1 older comment
Chris Horne on 31 Mar 2022
Yes, dx/dt = (sqrt(t+1)*cos(t^2)) / 2*(t+1) ) - 2*t*(sqrt(t+1)*sin(t^2))
with exact solution x = sqrt(1+t)*cos(t^2)

Torsten on 1 Apr 2022
% 2D SYSTEM SIN-COS EQUATIONS USING FORWARD EULER METHOD
% solve u'(t) = F(t,u(t)) where u(t)= sqrt(t+1)*cos^2(t),sqrt(t+1)*sin^2(t) ;
% % Initial conditions and setup
% x(t)= sqrt(t+1)*cos(t).^2
% y(t)= sqrt(t+1)*sin(t).^2
t(1)=0; % initial t
dt=0.005; % time step
nn=1000; % number of time steps in radians
t = zeros(1,nn);
x = zeros(1,nn);
y = zeros(1,nn);
x(1) = 1;
y(1) = 0; % initial x,y
for i=1:nn-1 % loop counter
dx=sqrt(t(i)+1)*cos(t(i)^2)./(2*(t(i)+1)) - 2*t(i)*(sqrt(t(i)+1)*sin(t(i)^2));%x' equation
dy=sqrt(t(i)+1)*sin(t(i)^2)./(2*(t(i)+1)) + 2*t(i)*(sqrt(t(i)+1)*cos(t(i)^2)); %y' equation
x(i+1)=x(i)+dt*dx; % Find new x
y(i+1)=y(i)+dt*dy; % Find new y
t(i+1)=t(i)+dt; % Find new t
end
% x Exact value without Euler
tt = linspace(t(1),t(end),numel(t));%set the tt,x vectors to same length
x_exact = sqrt(1+tt).*cos(tt.^2);
% Plot x vs. t with exact solution
figure
%t = linspace(0,1000,length(x));%set the t,x vectors to same length
plot(t,x,'b') % plot approximation in blue
hold on
title('2D System sqrt(1+t)cos equation vs exact soln')
xlabel('t')
ylabel('x')
plot(tt,x_exact,'r') % plot exact in red
legend('x','exact')
hold off;
figure
% y Exact value without Euler
tt = linspace(t(1),t(end),numel(t));%set the tt,y vectors to same length
y_exact = sqrt(1+tt).*sin(tt.^2);
% Plot y vs. t with exact solution
plot(t,y,'b')
hold on
title('2D System sqrt(1+t)sin equation vs exact soln')
xlabel('t')
ylabel('y')
%t = linspace(0,1000,length(y_exact));%set the t,y vectors to same length
plot(tt,y_exact,'r')
hold on
legend('y','exact')
hold off
Chris Horne on 1 Apr 2022
Torsten, Thanks!

James Tursa on 31 Mar 2022
Edited: James Tursa on 31 Mar 2022
You've only got one scalar differential equation, so I don't understand why you think you need two variables x and y along with independent variable t to handle this. You only need the x and t that are present in your original equation. E.g., something like this outline for the Euler part based on your posted code:
% % Initial conditions and setup
nn = 10000; % number of time steps in radians
t = zeros(1,nn);
x = zeros(1,nn);
t(1) = 0; % initial t
x(1) = 1; % initial x
dt = 0.005; % time step
for i=1:nn-1 % loop counter
dx = sqrt(t(i)+1)*cos(t(i)^2)/(2*(t(i)+1)) - 2*t(i)*(sqrt(t(i)+1)*sin(t(i)^2)); %x' equation
x(i+1) = x(i) + dt*dx; % Find new x
t(i+1) = t(i) + dt; % Find new t
end
So, only one derivative equation for the x. And all of those i's get replaced with t(i)'s. Also I have added a pair of parentheses to ensure that the entire 2*(t(i)+1) is in the denominator where it belongs.
For the exact solution, again you need only the one scalar x equation. And you need to make sure you are using the same t scale. E.g.,
tt = linspace(t(1),t(end),numel(t));
Finally, it is not clear to me what your intended final time is. In your Euler code it seems to be 9999*dt, but in your exact solution code it seems to be 10000. You need to make sure these are the same in order for your plots to match up properly.
Chris Horne on 31 Mar 2022
James, Thanks for the tip on the t(i). Then I did not expect the plot results to be so much different. The 'cos' code is same as the 'sin' code, more or less. The exact solution for y is much different than x.
% 2D SYSTEM SIN-COS EQUATIONS USING FORWARD EULER METHOD
% solve u'(t) = F(t,u(t)) where u(t)= sqrt(t+1)*cos^2(t),sqrt(t+1)*sin^2(t) ;
% % Initial conditions and setup
% x(t)= sqrt(t+1)*cos(t).^2
% y(t)= sqrt(t+1)*sin(t).^2
t(1)=0; % initial t
x(1)=1; y(1)=0; % initial x,y
dt=0.005; % time step
nn=1000; % number of time steps in radians
t = zeros(1,nn);
x = zeros(1,nn);
for i=1:nn-1 % loop counter
dx=sqrt(t(i)+1)*cos(t(i)^2)./(2*(t(i)+1)) - 2*t(i)*(sqrt(t(i)+1)*sin(t(i)^2));%x' equation
dy=sqrt(t(i)+1)*sin(t(i)^2)./(2*(t(i)+1)) - 2*t(i)*(sqrt(t(i)+1)*cos(t(i)^2)); %y' equation
x(i+1)=x(i)+dt*dx; % Find new x
y(i+1)=y(i)+dt*dy; % Find new y
t(i+1)=t(i)+dt; % Find new t
end
% x Exact value without Euler
tt = linspace(t(1),t(end),numel(t));%set the tt,x vectors to same length
x_exact = sqrt(1+tt).*cos(tt.^2);
% Plot x vs. t with exact solution
figure
t = linspace(0,1000,length(x));%set the t,x vectors to same length
plot(t,x,'b') % plot approximation in blue
hold on
title('2D System sqrt(1+t)cos equation vs exact soln')
xlabel('t')
ylabel('x')
plot(t,x_exact,'r') % plot exact in red
legend('x','exact')
hold off;
figure
% y Exact value without Euler
tt = linspace(t(1),t(end),numel(t));%set the tt,y vectors to same length
y_exact = sqrt(1+tt).*sin(tt.^2);
% Plot y vs. t with exact solution
plot(t,y,'b')
hold on
title('2D System sqrt(1+t)sin equation vs exact soln')
xlabel('t')
ylabel('y')
t = linspace(0,1000,length(y_exact));%set the t,y vectors to same length
plot(t,y_exact,'r')
hold on
legend('y','exact')
hold off;

R2021b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!