向前欧拉法求第一阶微分方程在Matlab

Question

我写是应该解决的IVP DU/DX = -5000（U（T） - cos（t））的 - 这Matlab程序SIN（t）的其中u（ 0）= 1。我的确切解决方案应该是u（t）= cos（t），但是我从我的代码中得到的欧拉向前的解决方案与我应该做的和我计算的结果相比是巨大的，但我不确定我已经走到哪里了在我的代码错误。你能找到我的错误吗？

function main 

dt=5; 
u0 = 1; 
n=50; 

[T,U] = euler(dt, u0, n); 

uexact = cos(T); 

plot(T,U) 

hold on 

plot(T, uexact, 'r') 

end 

function [T,U]= euler(dt, u0, n) 

R= dt/n; 
T=zeros(1,n+1); 
U=zeros(1,n+1); 
U(1)=u0; 
T(1) = 0; 


for j=1:n 

U(j+1)= U(j)+ R*(rhs(T(j),U(j))); 
T(j+1)= T(j) + R; 
end 

end 


function dP = rhs(t, P) 
P = zeros(1,1); 

dP = (-5000)*(P - cos(t)) - sin(t); 

end 

Answer 1

您不必在函数中设置P =零（1,1），该函数近似于公式的导数。

此外，你所面临的问题是[正向欧拉法的数值不稳定性] [1]。你需要一个非常小的时间步骤来使它收敛（因为在函数dP内部有很大的P系数）。

function main 

dt=5; 
u0 = 1; 
n=50000; % increase the number or subintervals (smaller time step) to get stability 

[T,U] = euler(dt, u0, n); 

uexact = cos(T); 

plot(T,U,'bs') 

hold on 

plot(T, uexact, 'r') 

end 

function [T,U]= euler(dt, u0, n) 

R= dt/n; 
T=zeros(1,n+1); 
U=zeros(1,n+1); 
U(1)=u0; 
T(1) = 0; 

for j=1:n 
% Implicit method converges 
% U(j+1) = (U(j) - R*(-5000)*cos(T(j)) - R*sin(T(j)))/(1 - R*(-5000)); 
U(j+1)= U(j)+ R*(rhs(T(j),U(j))); 
T(j+1)= T(j) + R; 
end 

end 


function dP = rhs(t, P) 
%P = zeros(1,1); %% It must not be here 

dP = (-5000)*(P - cos(t)) - sin(t); 

end 


    [1]: http://en.wikipedia.org/wiki/Euler_method