I need to solve a cauchy problem with the Euler Method.

24 次查看（过去 30 天）
Stefano Baioli
Stefano Baioli 2020-5-25
评论： Ameer Hamza 2020-5-25
Hi all. I was asked to solve this problem by my teacher:
I have to write a function that solves this cauchy problem with the Eulerian method, using an h (step size) of 0.25, in the interval [0,2].
I am then told:
"knowing that the exact solution of the equation is y(x)=cos(x)+sin(x), calculate the error committed in estimating y(2)."
I have to also make a graph of both the true solution (in one colour) and the solution made with the eulerian method in another colour.
What i've done so far (Before reaching the limit of my capabilities!):
I've written this script:
% Eulerian Method
% Setup and initial conditions
h = 0.25; % step size
x = 0:h:5; % interval of x
y = zeros(size(x)); % allocate the resulting y
y(1) = 1; % inizial value of y
n = numel(y); % the number of y values
% The loop that solves the differential equation:
for i=1:n-1
f = (y(i)-2*sin(x(i)))
y(i+1) = y(i) + h * f;
end
And i've managed to plot it into a graph but i have no idea of how to complete the tasks properly.
Please Help!

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

采纳的回答

Ameer Hamza
Ameer Hamza 2020-5-25
Following shows how to find true solution using symbolic toolbox and plot both of them on same graph
% Eulerian Method
% Setup and initial conditions
h = 0.25; % step size
x = 0:h:5; % interval of x
y = zeros(size(x)); % allocate the resulting y
y(1) = 1; % inizial value of y
n = numel(y); % the number of y values
% The loop that solves the differential equation:
for i=1:n-1
f = (y(i)-2*sin(x(i)));
y(i+1) = y(i) + h * f;
end
syms Y(X)
eq = diff(Y) == Y - 2*sin(X);
ic = Y(0) == 1;
sol = dsolve(eq, ic);
y_sol = double(subs(sol, X, x));
hold on
plot(x, y_sol, 'r', 'DisplayName', 'true solution');
plot(x, y, 'b', 'DisplayName', 'Euler method');
legend
  5 个评论
显示 3更早的评论
John D'Errico
John D'Errico 2020-5-25
编辑：John D'Errico 2020-5-25
I would suggest it is a bad idea to use interp1 there, certainly so if you use the default method for interp1. The default is LINEAR interpolation. In that case, on a coarse interval, a significant fraction of the error you might report could just be that of interp1!
In context of what you are doing, it just so happens that x==2 just happens to be one of the points where you predicted the solution. In that case, linear interpolation is meaningless, since you would just want to see the value at the point of interest, and interp1 will do exactly the same. However, suppose the question had been to predict the error at some intermediate point? In that case curvature of the result would be a bad thing, since the connect-the-dots linear interpolation done inside interp1 will hurt you.
In the specific case here where you just need the error at x==2, you could have done just this:
Error = y(9) - (sin(x) + cos(x));
In a more general case, perhaps where you were asked to provide the error at say x=2.13, thus intermediate to your data, now you might more accurately have done:
xerr = 2.13;
ytrue = @(x) sin(x) + cos(x);
Error = interp1(x,y,xerr,'spline') - ytrue(xerr);
At least here the error of interpolation will be reduced due to the use of a spline, and only for the vector y. See here that I used the true function to predict the value at the point of interest, not an interpolation at all.
Always try to avoid adding additional sources of error to what you do.
An even better choice would have been to carefully write way-points into your code, thus making sure the ODE solver predicts a value at exactly some points of interst. This requires more carefully written code of course.
Ameer Hamza
Ameer Hamza 2020-5-25
Yes, interp1 will introduce its own error, however, in the previous case, I ignored it because the query point xq=2 also lies in vector x. However, in the general case, it can fail. Following code is for multiple values of h and calculates error using the method suggested by John
% Eulerian Method
% Setup and initial conditions
H = [0.25 0.1 0.05 0.01 0.005]; % step size
Y_euler = cell(1,numel(H));
for k=1:numel(H)
h = H(k);
x = 0:h:5; % interval of x
y = zeros(1, numel(x)); % allocate the resulting y
y(1) = 1; % inizial value of y
n = numel(y); % the number of y values
% The loop that solves the differential equation:
for i=1:n-1
f = (y(i)-2*sin(x(i)));
y(i+1) = y(i) + h * f;
end
Y_euler{k} = [x; y].';
end
syms Y(X)
eq = diff(Y) == Y - 2*sin(X);
ic = Y(0) == 1;
sol = dsolve(eq, ic);
y_sol = double(subs(sol, X, x));
hold on
plot(x, y_sol, 'r', 'DisplayName', 'true solution');
for k=1:numel(H)
plot(Y_euler{k}(:,1), Y_euler{k}(:,2), 'b--', 'DisplayName', ['Euler method h=' num2str(H(k), '%.3f')]);
end
legend
err = zeros(1,numel(H));
for k=1:numel(H)
err(k) = interp1(Y_euler{k}(:,1), Y_euler{k}(:,2), 2) - subs(sol, X, 2);
end

请先登录，再进行评论。

更多回答（0 个）

类别

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

标签

产品


版本

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by