ODE solution using Fourier series

7 visualizzazioni (ultimi 30 giorni)
Muhammad Usman
Muhammad Usman il 16 Mag 2021
Risposto: Jan il 16 Mag 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 Commento
Jan
Jan il 16 Mag 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.

Accedi per commentare.

Risposte (1)

Jan
Jan il 16 Mag 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.

Community Treasure Hunt

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

Start Hunting!

Translated by