Sine equation in Euler Method ODE

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')

[James Tursa]

Can you clarify what the differential equation is you are solving? Maybe post an image of the equation? There's lots of things obviously wrong with your posted code, but before getting into all of that it would be best to verify the differential equations and initial conditions we are working with first.

[Chris Horne]

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)

Thanks for your interest

Answer 1

[Torsten]

% 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]

Thanks, My plots are way too different. Surprising since dx/dt and dy/dt are the same except for cos in x' and sin in y'.

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 diffeq

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 diffeq

[Torsten]

You made some mistakes in your code.

I corrected them, and now the plots for approximate and analytical solution coincide.

[Chris Horne]

Torsten, Thanks!

Answer 2

[James Tursa]

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]

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;