Sine equation in Euler Method ODE

8 visualizaciones (últimos 30 días)
Chris Horne
Chris Horne el 31 de Mzo. de 2022
Comentada: Chris Horne el 1 de Abr. de 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 comentarios
James Tursa
James Tursa el 31 de Mzo. de 2022
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
Chris Horne el 31 de Mzo. de 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)
Thanks for your interest

Iniciar sesión para comentar.

Respuesta aceptada

Torsten
Torsten el 1 de Abr. de 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
  3 comentarios
Torsten
Torsten el 1 de Abr. de 2022
Editada: Torsten el 1 de Abr. de 2022
You made some mistakes in your code.
I corrected them, and now the plots for approximate and analytical solution coincide.
Chris Horne
Chris Horne el 1 de Abr. de 2022
Torsten, Thanks!

Iniciar sesión para comentar.

Más respuestas (1)

James Tursa
James Tursa el 31 de Mzo. de 2022
Editada: James Tursa el 31 de Mzo. de 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.
  1 comentario
Chris Horne
Chris Horne el 31 de Mzo. de 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;

Iniciar sesión para comentar.

Categorías

Más información sobre Programming en Help Center y File Exchange.

Productos


Versión

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by