Модифицированный метод Эйлера.
Для определения 
воспользуемся обычной формулой Эйлера, а для определения 
будет естественным использовать уточненную формулу, т.е.
(Л4.3)
Пример соответствующей M-функции (ниже задано 
, 
, 
):
function consol
s=input('s=');g=input('g=');L=input('L=');n=input('n=');
h=L/n;c=0.02*(s+g);ej=1/c;
x=0:h:L;z=zeros(size(x));y=z;
for i=1:n
z(i+1)=z(i)+h*f(x(i),z(i));
y(i+1)=y(i)+h*(z(i)+z(i+1))/2;
end
ixyz=[(0:n);x;y;z];
fprintf('\n'),fprintf('%3s %4s %6s %6s','i','x','y','z')
fprintf('\n'),fprintf('%3d %4.1f %6.4f %6.4f\n',ixyz)
plot(x,y,'.-',x,z,'*-'),grid on
legend('y=y(x)','z=y''(x)',0),title('modificyrovannyj metod Eilera')
function dzdx=f(x,z)
dzdx=m(x)/ej*sqrt((1+z^2)^3);
end
function M=m(x)
M=1/sqrt((1+(c*x)^2)^3);
end
end
Результаты счета:
s=12
g=3
L=1
n=10
i x y z
0 0.0 0.0000 0.0000
1 0.1 0.0015 0.0300
2 0.2 0.0060 0.0600
3 0.3 0.0135 0.0900
4 0.4 0.0240 0.1200
5 0.5 0.0375 0.1500
6 0.6 0.0540 0.1800
7 0.7 0.0735 0.2100
8 0.8 0.0960 0.2400
9 0.9 0.1215 0.2700
10 1.0 0.1500 0.3000
>>
2. Использование встроенной функции ode45, предназначенной для решения задачи Коши для систем обыкновенных дифференциальных уравнений 1-го порядка вида 
.
В данной функции реализуются одношаговые методы Рунге-Кутта 4-го и 5-го порядков. Это классический подход, рекомендуемый для начальной пробы решения. Во многих случаях (если система уравнений не является жесткой) можно получить удовлетворительные результаты. Таким образом, данную функцию применяют довольно часто, особенно тогда, когда неизвестны характеристики задачи.
Предварительно представим систему (Л4.2) в виде
, (Л4.3)
где 
, 
, 
,
, 
.
При обращении к ode45 входными параметрами будут дескриптор имени функции, реализующей вычисление вектор-функции 
, массив, содержащий значения левого и правого конца диапазона изменения переменной 
, а также вектор начальных значений 
. Выходными параметрами будут массивы 
и 
, где
, 
, 
, 
.
Пример соответствующей M-функции (ниже задано , ):
function consol_ode45
s=input('s=');
g=input('g=');
c=0.02*(s+g);
EJ=1/c;
x0xL=input('[x0 xL]=');
z0=input('[z10 z20]=');
[X Z]=ode45(@F,x0xL,z0);
n=length(X);
fprintf('\n'),fprintf('%3s %6s %6s %6s','i','x','y','y'''),fprintf('\n')
for i=1:n
fprintf('%3d %6.4f %6.4f %6.4f\n',i,X(i),Z(i,:))
end
plot(X,Z(:,1),'.-',X,Z(:,2),'*-'),grid on
legend('y=y(x)','z=y''(x)',0),title('ode45')
function dz=F(x,z)
dz=zeros(2,1);
dz(1)=z(2);
dz(2)=M(x)/EJ*sqrt((1+z(2)^2)^3);
end
function m=M(x)
m=1/sqrt((1+(c*x)^2)^3);
end
end
Результаты счета:
s=12
g=3
[x0 xL]=[0 1]
[y0 z0]=[0 0]
i x y y'
1 0.0000 0.0000 0.0000
2 0.0002 0.0000 0.0001
3 0.0003 0.0000 0.0001
4 0.0005 0.0000 0.0002
5 0.0007 0.0000 0.0002
6 0.0015 0.0000 0.0005
7 0.0023 0.0000 0.0007
8 0.0032 0.0000 0.0010
9 0.0040 0.0000 0.0012
10 0.0082 0.0000 0.0025
11 0.0124 0.0000 0.0037
12 0.0166 0.0000 0.0050
13 0.0208 0.0001 0.0062
14 0.0417 0.0003 0.0125
15 0.0626 0.0006 0.0188
16 0.0836 0.0010 0.0251
17 0.1045 0.0016 0.0313
18 0.1295 0.0025 0.0388
19 0.1545 0.0036 0.0463
20 0.1795 0.0048 0.0538
21 0.2045 0.0063 0.0613
22 0.2295 0.0079 0.0688
23 0.2545 0.0097 0.0763
24 0.2795 0.0117 0.0838
25 0.3045 0.0139 0.0913
26 0.3295 0.0163 0.0988
27 0.3545 0.0188 0.1063
28 0.3795 0.0216 0.1138
29 0.4045 0.0245 0.1213
30 0.4295 0.0277 0.1288
31 0.4545 0.0310 0.1363
32 0.4795 0.0345 0.1438
33 0.5045 0.0382 0.1513
34 0.5295 0.0421 0.1588
35 0.5545 0.0461 0.1663
36 0.5795 0.0504 0.1738
37 0.6045 0.0548 0.1813
38 0.6295 0.0594 0.1888
39 0.6545 0.0643 0.1963
40 0.6795 0.0693 0.2038
41 0.7045 0.0744 0.2113
42 0.7295 0.0798 0.2188
43 0.7545 0.0854 0.2263
44 0.7795 0.0911 0.2338
45 0.8045 0.0971 0.2413
46 0.8295 0.1032 0.2488
47 0.8545 0.1095 0.2563
48 0.8795 0.1160 0.2638
49 0.9045 0.1227 0.2713
50 0.9284 0.1293 0.2785
51 0.9522 0.1360 0.2857
52 0.9761 0.1429 0.2928
53 1.0000 0.1500 0.3000
>>
>> 
3. Использование встроенной функции ode15s, предназначенной для решения задачи Коши для систем обыкновенных дифференциальных уравнений 1-го порядка вида .
В данной функции реализуется многошаговый метод переменного порядка (от одного до пяти), использующий формулы численного дифференцирования «назад». Это адаптивный метод, рекомендуемый к применению в случае, когда ode45 не обеспечивает решения (система уравнений жесткая). При этом, входные и выходные параметры при обращении к ode15s полностью совпадают с соответствующими параметрами функции ode45 .
Пример соответствующей M-функции (ниже задано , ):
function consol_ode15s
s=input('s=');
g=input('g=');
c=0.02*(s+g);
EJ=1/c;
x0xL=input('[x0 xL]=');
z0=input('[z10 z20]=');
[X Z]=ode15s(@F,x0xL,z0);
n=length(X);
fprintf('\n'),fprintf('%3s %6s %6s %6s','i','x','y','y'''),fprintf('\n')
for i=1:n
fprintf('%3d %6.4f %6.4f %6.4f\n',i,X(i),Z(i,:))
end
plot(X,Z(:,1),'.-',X,Z(:,2),'*-'),grid on
legend('y=y(x)','z=y''(x)',0),title('ode15s')
function dz=F(x,z)
dz=zeros(2,1);
dz(1)=z(2);
dz(2)=M(x)/EJ*sqrt((1+z(2)^2)^3);
end
function m=M(x)
m=1/sqrt((1+(c*x)^2)^3);
end
end
Результаты расчета:
s=12
g=3
[x0 xL]=[0 1]
[y0 z0]=[0 0]
i x y y'
1 0.0000 0.0000 0.0000
2 0.0021 0.0000 0.0006
3 0.0041 0.0000 0.0012
4 0.0062 0.0000 0.0019
5 0.0152 0.0000 0.0046
6 0.0242 0.0001 0.0073
7 0.0332 0.0002 0.0100
8 0.0422 0.0003 0.0127
9 0.0690 0.0007 0.0207
10 0.0957 0.0014 0.0287
11 0.1225 0.0023 0.0367
12 0.1492 0.0033 0.0448
13 0.2492 0.0093 0.0748
14 0.3492 0.0183 0.1048
15 0.4492 0.0303 0.1348
16 0.5492 0.0453 0.1648
17 0.6492 0.0632 0.1948
18 0.7492 0.0842 0.2248
19 0.8492 0.1082 0.2548
20 0.9492 0.1352 0.2848
21 1.0000 0.1500 0.3000
>>
