Trying to change the output vector for Heun's non-self-starting

1 次查看（过去 30 天）
Paul Toepp
Paul Toepp 2022-4-13
回答： Shivam 2023-10-16
I have created a program to run for non-self-starting Heun's method for a given problem in class. But what I am noticing is that my "y" output has one extra 0 at the start that id like to remove.
This is the given problem
And this is my code
f= @(t,y)(sin(y)/y+2*t);
a = input('Enter LOWER limit of domain, a= ');
b = input('Enter UPPER limit of domain, b= ');
h = input('Enter step size, h= ');
t=a:h:b;
y=zeros(size(t));
y(2) = input('Enter Initial condition of y(0),= ');
y(1) = input('Enter Initial condition of y(i-1),= ');
n=numel(t);
for i=2:n
nil(i-1)=y(i-1)+(2*h*f(t(i-1),y(i)));
y(i+1)=y(i)+(h/2)*(f(t(i),nil(i-1))+f(t(i-1),y(i)));
end
Output is as follows (im only interested in the "t" andn "y") Now I very well could be wrong in this but given the initial condition of the problem, should my first "y" output be 1 and not 0? and if so, how can i fix that.
  1 个评论
Johan
Johan 2022-4-13
The first two values of y are your initial conditions:
y(2) = input('Enter Initial condition of y(0),= ');
y(1) = input('Enter Initial condition of y(i-1),= ');
According to your problem they should be 0 and 1 if I understand correctly.
I'm not familiar with non self starting Heun's method but the solution may be to shift your vector t by -1 step ?
t=(a-h):h:b;

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Shivam
Shivam 2023-10-16
Hi Paul,
As per my understanding, you are having an issue filling output vector 'y' while solving given differential equation with non-self-starting Heun's method.
You can follow the below workaround to resolve the issue:
  • Modify vector t to have its first value as -1.
t = [-1, a:h:b];
  • Modify for loop to calculate the correct predictor and corrector following the non-self-starting Heun's method.
% ..
% ...
n=numel(t);
for i=2:n-1
nil = y(i-1)+f(t(i), y(i))*2*h; % Predictor
y(i+1) = y(i)+(h/2)*(f(t(i), y(i))+f(t(i+1), nil)); % Corrector
end
The first value of vector y corresponds to y(t(1)) ~= y(-1) = 0.
In case you need the output vector y with no output value corresponding to t = -1, you can modify y after above calculations:
y = y(:, 2:end)
I hope it helps in resolving the issue.

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by