Sine equation in Euler Method ODE

8 visualizzazioni (ultimi 30 giorni)
Chris Horne
Chris Horne il 31 Mar 2022
Commentato: Chris Horne il 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 Commenti
James Tursa
James Tursa il 31 Mar 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 il 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)
Thanks for your interest

Accedi per commentare.

Risposta accettata

Torsten
Torsten il 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
  3 Commenti
Torsten
Torsten il 1 Apr 2022
Modificato: Torsten il 1 Apr 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 il 1 Apr 2022
Torsten, Thanks!

Accedi per commentare.

Più risposte (1)

James Tursa
James Tursa il 31 Mar 2022
Modificato: James Tursa il 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.
  1 Commento
Chris Horne
Chris Horne il 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;

Accedi per commentare.

Categorie

Scopri di più su Programming in Help Center e File Exchange

Prodotti


Release

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by