Euler's Method for stiff ODE
13 ビュー (過去 30 日間)
古いコメントを表示
Having problem with my code, I'm trying to make my graph (the middle) look like the bottom graph. I feel like it is something with the initial inputs, I pretty sure my for loop is correct. Thanks for the help I appreciate it!
if true
% %%Explicit Euler
h = 5;
x = [0:h:20];
yb = zeros(size(x));
yb(1) = 1;
yb2(1) = 1;
mu = 1;
for n = 1:length(x)-1
yb(n+1) = yb(n) + h*yb2(n);
yb2(n+1) = yb2(n) + h*(mu*(1-yb(n)^2)*yb2(n) - yb(n));
end
%%Plot
figure(1); clf(1);
plot(x,yb,'k');
hold on
plot(x,yb2,'or');
title(['Van der Pol Equation for \mu = ',num2str(mu), 'with Euler']);
xlabel('Time');
ylabel('y(t)');
legend('y_{1}', 'y_{2}','location','northeast')
end
4 件のコメント
採用された回答
Jim Riggs
2018 年 3 月 1 日
編集済み: Jim Riggs
2018 年 3 月 1 日
Torsten is right. There is no way you will get a solution with a large step size (h=5). I checked out your code and it works fine for mu=1 with a step size of h=0.1
I would recommend the following two changes;
1) add the command
clear x yb yb2
at the start of the script.
2) I would also recommending pre-allocating the yb2 vector. you can do this by changing
yb2(1) = 1
to
yb2 = yb.
this will make yb2 and yb the same size, and they will both have a 1 in the first position.
Now, as you increase the value of mu, the system becomes more stiff. For example, using a step size of h=0.1 I can get a solution for mu =2 and mu=3, but at mu=4 the function blows up after about 14 seconds. (it is numerically unstable). Now you have to reduce the step size to improve the calculation stability. In order to get a solution with mu=1000, and to execute for 3000 seconds without blowing up, the step size must be extremely small (10^-6 or even smaller). My problem is that this is taxing the capabilities of my computer and I cannot find the solution.
A major part of the problem is the use of the Euler integration method. This method is not well suited for stiff systems (perhaps that is what this lesson is about). For this problem, I would switch to a Runge-Kutta integrator and try my favorite [5,6] Runge Kutta method. This would allow the use of a much larger step size (probably 3 orders of magnitude) compared to what is required for Euler integration.
2 件のコメント
Jim Riggs
2018 年 3 月 2 日
編集済み: Jim Riggs
2018 年 3 月 2 日
This will work for mu=1, since it results in a step size of 0.04. Bit it will not work for significantly larger values of mu. As mu increases, you still need to decrease the step size. You will have to experiment to find acceptable values. And when you increase the timespan, the step size gets bigger. This is also not desirable.
その他の回答 (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!