Chebyshev Differentiation Matrix to solve ODE
显示 更早的评论
Trying to solve dy/dx = 2xy, y(1) = 1, using Chebyshev differentiation matrix
The exact solution is y = exp(x.^2 -1);
Heres what I have:
% compute the chebyshev differentiation matrix and x-grid
[D1,x] = cheb(N);
% boundary condition (I don't know if this is correct?)
D1 = D1(1:end-1,1:end-1); x = x(1:end-1);
% compute the derivatives at x (i.e. at the chebyshev grid points)
f = 2*x.*ones(size(x));
% solve
u = D1\f;
% set the boundary condition
u = [u;1];
Where cheb.m is from Trefethen (spectral methods in matlab)
function [D,x] = cheb(N)
% check the base case
if N == 0; D = 0; x = 1; return; end
% create the Chebyshev grid
x = cos(pi*(0:N)/N)';
c = [2; ones(N-1,1);2].*(-1).^(0:N)';
X = repmat(x,1,N+1);
dX = X-X';
D = (c*(1./c)')./(dX+(eye(N+1)));
D = D - diag(sum(D'));
This solution (u = D1\f) does not match the exact solution at all.
I think what I have is close ... Any help would be awesome. Thanks in advance!
采纳的回答
Although this question is more than a decade old, there still seems to be some interest in this. I was suprised to see that all the answers use the solution y = exp(x.^2 -1) which is what we are supposed to find by solving the ODE dy/dx = 2xy, y(1) = 1 numerically.
clear
clc
% Solve: dy/dx = 2xy, y(x), y(1) = 1
N = 20;
% compute the chebyshev differentiation matrix and x-grid
[D,x] = cheb(N); %dom=[-1,1]
% Domain is important as we will be solving for y only inside this domain.
% The ODE problem can now be written as the following linear problem with
% one constraint. This constraint arises due to the boundary condition.
% D*y_est == 2*x.*y_est
% y(x=1) = 1.
y_est = nan(N+1,1);
y_est(end) = 1; % boundary condition
% Note that we have N unknowns (values of y at x(1:N)) but the system above has N+1 equations.
% To solve this we utilize the boundary condition (constraint) to reduce
% the system to the form q = Ap with size(A) = [N,N] (N equations).
% Solving this using mldivide (\ operator) gives us
y_est(1:N) = (D(1:N,1:N)-2*diag(x(1:N)))\(-D(1:N,N+1)*1);
% Remark that we didn't use y = exp(x.^2 -1). In the figure below we can
% see that the solution to the ODE we calculated matches the exact solution
% y = exp(x.^2 -1) well at the x-grid points.
% Plotting
plot(x,y_est,'s',x,exp(x.^2-1))
xlabel("x")
ylabel("y, y_{est}")
legend('Chebyshev','Exact')
%%
% Checking convergence. Here we can see how the numerically computed
% solution improves as we increase number of x-grid points N.
figure
for N = 2:30
% compute the chebyshev differentiation matrix and x-grid
[D,x] = cheb(N); %dom=[-1,1]
% Domain is important as we will be solving for y only inside this domain.
% The ODE problem can now be written as the following linear problem with
% one constraint. This constraint arises due to the boundary condition
% y(x=1) = 1.
% D*y_est == 2*x.*y_est
y_est = nan(N+1,1);
y_est(end) = 1; % boundary condition
% Note that we have N unknowns (values of y at x(1:N)) but the system above has N+1 equations.
% To solve this we utilize the boundary condition (constraint) to reduce
% the system to the form q = Ap with size(A) = [N,N].
% Solving this using mldivide (\ operator) gives us
y_est(1:N) = (D(1:N,1:N)-2*diag(x(1:N)))\(-D(1:N,N+1)*1);
semilogy(N,norm(y_est-exp(x.^2-1)),'bo')
xlabel('N')
ylabel("||y-y_{est}||")
hold on
drawnow
end
function [D,x] = cheb(N)
% check the base case
if N == 0; D = 0; x = 1; return; end
% create the Chebyshev grid
x = cos(pi*(0:N)/N)';
c = [2; ones(N-1,1);2].*(-1).^(0:N)';
X = repmat(x,1,N+1);
dX = X-X';
D = (c*(1./c)')./(dX+(eye(N+1)));
D = D - diag(sum(D'));
end
1 个评论
John D'Errico
2024-4-4
编辑：John D'Errico
2024-4-4
I had to laugh myself when I was reading through some of the other answers, since I was observing the same thing. People seem to be trying to solve an ode using the known exact solution. The good news is, I could accept your answer.
更多回答（1 个）
The original post was right that the boundary condition was not set correct. Here's a right way to set the boundary condition:
% compute the chebyshev differentiation matrix and x-grid
[D,x] = cheb(N);
% compute the derivatives at x (i.e. at the chebyshev grid points)
f = 2*x.*exp(x.^2 - 1);
% boundary condition (correct way)
D(end,:) = 0;
D(end,end) = 1;
f(end) = 1;
% solve
u = D\f;
2 个评论
Lukgaf
2018-3-19
I used the cod but not calable. kindly help
Walter Roberson
2018-3-19
What value are you passing? What error message are you receiving?
You opened a Question on this topic, so the discussion should probably be there.
类别
在 帮助中心 和 File Exchange 中查找有关 Mathematics 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
