blob: ce01ffb75899d20f9a781772425709a871a7d04c (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
|
//clc()
//f'(x,y) = -2*x^3 + 12*x^2 -20*x + 8.5
//f(x,y) = -x^4 / 2 + 4*x^3 - 10*x^2 + 8.5*x + 1
h = 0.5;
x = 0:h:4;
y1 = -x^4 / 2 + 4*x^3 - 10*x^2 + 8.5*x + 1;
y(1) = 1;
for i=1:8
k1(i) = -2*x(i)^3 + 12*x(i)^2 -20*x(i) + 8.5;
x1(i) = x(i) + h/2;
k2(i) = -2*x1(i)^3 + 12*x1(i)^2 -20*x1(i) + 8.5;
k3(i) = -2*x1(i)^3 + 12*x1(i)^2 -20*x1(i) + 8.5;
x2(i) = x(i) + h;
k4(i) = -2*x2(i)^3 + 12*x2(i)^2 -20*x2(i) + 8.5;
y(i+1) = y(i) + (k1(i)+2*k2(i)+2*k3(i)+k4(i))*h/6;
e(i) = (y1(i) - y(i))*100/y1(i);
end
disp("f(x,y) = -2*x^3 + 12*x^2 -20*x + 8.5")
disp(y(1:9),"y by fourth order Ralston RK")
disp("f(x,y) = 4*exp(0.8*x) - 0.5*y")
x = 0:h:0.5;
y(1) = 2;
k1 = 4*(exp(0.8*x(1)))-0.5*y(1);
x1 = x(1) + 0.5*h;
y1 = y(1) + 0.5*k1*h;
k2 = 4*exp(0.8*x1) - 0.5*y1;
y2 = y(1) + 0.5*k2*h;
k3 = 4*exp(0.8*x1) - 0.5*y2;
x1 = x(1) + h;
y1 = y(1) + k3*h;
k4 = 4*exp(0.8*x1) - 0.5*y1;
y(2) = y(1) + (k1+2*k2+2*k3+k4)*h/6;
disp(y(1:2),"y = ")
|
