数值分析第九章作业

第九章
欧拉法

function [X,Y] = Euler(a,b,h,f)
n = abs((a-b)/h);
X =a:h:b;
Y =zeros(1,n);
Y(1) = 1;

for i = 1:n
    Y(i+1) = Y(i) +h*  f(X(i),Y(i));
end

改进的欧拉法

function [X,Y] = Euler_gaijin(a,b,h,f)
n = abs((a-b)/h);
X =a:h:b;
Y =zeros(1,n);
Y(1) = 1;

for i = 1:n
    % 预测
    Y(i+1) = Y(i) +h*  f(X(i),Y(i));
    % 校正
    Y(i+1) = Y(i) + h/2 * (f(X(i),Y(i)) + f(X(i+1),Y(i+1)));
end

主程序

clc 
clear
%% 准备
a = 1 ;
b = 2;
h = 0.1;
f = @(x,y)(1/x^2)-y/x;

%% 欧拉法
[X1,Y1] = Euler(a,b,h,f);
subplot(4,3,1),plot(X1,Y1)
[X2,Y2] = Euler(a,b,0.025,f);
subplot(4,3,2),plot(X2,Y2)
[X3,Y3] = Euler(a,b,0.01,f);
subplot(4,3,3),plot(X3,Y3)
%% 改进的欧拉法
[X4,Y4] = Euler_gaijin(a,b,h,f);
subplot(4,3,4),plot(X4,Y4)
[X5,Y5] = Euler_gaijin(a,b,0.025,f);
subplot(4,3,5),plot(X5,Y5)
[X6,Y6] = Euler_gaijin(a,b,0.01,f);
subplot(4,3,6),plot(X6,Y6)
subplot(4,3,[7,12]),
plot(X1,Y1)
hold on
plot(X2,Y2)
hold on
plot(X3,Y3)
hold on
plot(X4,Y4)
hold on
plot(X5,Y5)
hold on
plot(X6,Y6)

结果