我正在嘗試為 Forward Euler 方法撰寫 MATLAB 代碼，但我認為輸出不完全正確。這是問題：
撰寫實作前向歐拉方案的 Matalb 函式來求解 p ODE 的泛型系統。您的函式必須符合以下介面函式 [tHist ,uHist] = ForwardEuler(f,tspan ,u0,nsteps) 并具有以下輸入：
? (1) 的函式f 的函式句柄f；
? 一個 2 向量 tspan，具有分量 t0 和 T；
?ap ×1 向量u0，包含初始條件y0；
? 要采取的步驟 nteps 的數量；
以及以下輸出：
?an (Nh 1) ×1 向量 tHist，存盤近似時間 tn
?an (Nh 1) ×p 矩陣uHist，其第n 行存盤近似值un ≈y(tn) ∈R^p。
這是我的代碼：
function [tHist ,uHist] = ForwardEuler(f,tspan ,u0,nsteps)
%Inputs:
% f: a function handle for the function f
% tspan: a 2-vector tspan, with components t0 and T, the initial and final components respectively;
% u0 : a p ×1 vector containing the initial condition y0;
% nsteps : he number of steps to be taken;
% Outputs:
% tHist: a (Nh 1) ×1 vector storing the approximation times tn
% uHist: an (Nh 1) ×p matrix uHist, whose nth row stores the approximation un ≈y(tn) ∈ R^p.
T = [tspan(1) :nsteps: tspan(2)];
u = zeros(1,nsteps);
u(n) = u0;
h = (tspan(2) - tspan(1))/nsteps;
for i= 1: nsteps
u(i 1) = u(i) h*f(T(i),u(i));
tHist = [u(i 1)*T(i)];
uHist = [u(i 1)];
end
end
uj5u.com熱心網友回復：
一個最小的直接實作可能看起來像
function [tHist ,uHist] = ForwardEuler(f,tspan ,u0,nsteps)
tHist = linspace(tspan(1), tspan(2), nsteps 1);
h = tHist(2)-tHist(1);
uHist = zeros(nsteps 1,length(u0));
uHist(1,:)=u0;
for i = 1:nsteps
uHist(i 1,:) = uHist(i,:) h*f(tHist(i),uHist(i,:));
end%for
end%function
用圓方程測驗例程
[T,U] = ForwardEuler(@(t,u)[u(2), -u(1)], [0,6.5], [1,0], 60);
plot(U(:,1),U(:,2));
這產生了預期的向外螺旋。
uj5u.com熱心網友回復：
核心演算法似乎沒問題，只是必須更正初始值。
function [T ,u] = ForwardEuler(f,tspan ,u0,nsteps)
u = zeros(1,nsteps 1);%corrected to nsteps 1
u(1) = u0;%corrected n to 1
h = (tspan(2) - tspan(1))/nsteps;
T = [tspan(1) :h: tspan(2)];%corrected nsteps to h
for i= 1: nsteps
u(i 1) = u(i) h*f(T(i),u(i));
end
end
tHist 和 uHist 仍未完成，但沒有太多作業要做。
另一種解決方案是為每一步將 T 和 u 向量存盤在矩陣中。
function [tHist ,uHist] = ForwardEulerWithHist(f,tspan ,u0,nsteps)
tHist=zeros(nsteps 1,nsteps);
uHist=zeros(nsteps 1,nsteps);
for j=1:nsteps
[T,u] = ForwardEuler(f,tspan,u0,j);
tHist(1:(j 1),j) = T;
uHist(1:(j 1),j) = u;
end
end
然后您可以比較迭代平滑到最佳迭代的程度。示例比較步驟 20 的迭代與步驟 10 的迭代：
[tHist,uHist]=ForwardEulerWithHist(@(x,y)x y.^2,[0,0.5],0,20);
plot(tHist(:,20),uHist(:,20)), hold on;plot(tHist(1:11,10),uHist(1:11,10))
轉載請註明出處，本文鏈接：https://www.uj5u.com/caozuo/401345.html
上一篇：如何在引導程式中居中對齊輸入組