ODE solution using Fourier series
3 次查看(过去 30 天)
显示 更早的评论
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/619583/image.png)
Fourier series:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/619588/image.png)
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 个评论
Jan
2021-5-16
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.
回答(1 个)
Jan
2021-5-16
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))
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
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)
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.
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!