ODE solution using Fourier series

4 visualizaciones (últimos 30 días)
Muhammad Usman
Muhammad Usman el 16 de Mayo de 2021
Respondida: Jan el 16 de Mayo de 2021
Hello,
I am trying to solve first order ode () using Fourier Series. I wrote the code, but I was doing something, please idetify my mistake.
Fourier series:
by putting into original ODE. Here's the code:
clear all;clc;close all;
syms x
N = 3;
for j=1:(2*N+1)
xp(j) = simplify(-sym(pi) + (2*sym(pi)*(j-1))./(2.*N+1)); % Equidistance points
end
for k = 1:2*N+1
for l =1:N
a(k,l) = cos(l*xp(k));
b(k,l) = sin(l*xp(k));
end
end
one_col = ones(1,2*N+1);
sumab = double([a';b']);
A = (1+cos(x)).*double([one_col;sumab])'
x = (A\ones(7,1))
a0 = x(1);
an = x(2:N+1);
bn = x(N+2:2*N+1);
dA = diff(A);
% ODE
A_ODE = dA + A;
yxp = A_ODE*x
plot(xp,yxp);
hold on;
%% ODE45
[a,b] = ode45(@(a,b) -(1+cos(a))*b+1, [0 30],0); %to distinguish variable I choose x = a and y = b
plot(a,b)
Thanks
  1 comentario
Jan
Jan el 16 de Mayo de 2021
Please mention, why you assume, that there is a mistake. You d have this important information already. Then it is useful to share it with the ones, who want to help you.

Iniciar sesión para comentar.

Respuestas (1)

Jan
Jan el 16 de Mayo de 2021
syms x
N = 3;
for j=1:(2*N+1)
xp(j) = simplify(-sym(pi) + (2*sym(pi)*(j-1))./(2.*N+1)); % Equidistance points
end
for k = 1:2*N+1
for l =1:N
a(k,l) = cos(l*xp(k));
b(k,l) = sin(l*xp(k));
end
end
one_col = ones(1,2*N+1);
sumab = double([a';b']);
A = (1+cos(x)).*double([one_col;sumab])'
A = 
x = (A\ones(7,1))
x = 
Are you sure that this is wanted?
a0 = x(1);
an = x(2:N+1);
bn = x(N+2:2*N+1);
dA = diff(A);
% ODE
A_ODE = dA + A;
yxp = A_ODE*x
yxp = 
plot(xp,yxp);
Error using plot
Data must be numeric, datetime, duration or an array convertible to double.
hold on;
%% ODE45
[a,b] = ode45(@(a,b) -(1+cos(a))*b+1, [0 30],0); %to distinguish variable I choose x = a and y = b
plot(a,b)
Your code does not contain meaningful comments. Guessing its intention is not reliable. Therefore it is hard to estimate, if something is wanted or a mistake.

Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by