matlab中四阶龙格库塔算法、欧拉算法和改进的欧拉算法的总结

学习了一天matlab，把自己从书上学到的知识写出来。

重要的概念

常微分方程 （ordinary differential equation，简称ODE）是未知函数中只含有一个自变量的微分方程。

偏微分方程（partial differential equation，简称PDE）是未知函数是多元函数的微分方程。

CAE求解方法一般有两种

显式（Explicit）第n步结果可以从n-1, n-2, ......1步的结果直接推导出来，迭代时每步的计算量很小，但迭代增量也有限制，不能太大，否则会出现发散。

隐式（Implicit）第n步的计算结果不能直接从前面的结果推导出来，必须做进一步的求解，这样，迭代时每步的计算量很大。

泰勒公式

若 f(x) 在 x = 0 点直到 n + 1 阶连续导数，则

称为 f(x) 在 x = 0 点关于 x 的幂函数展开式，又称为 Taylor 公式，式中Rn（x）叫做 Lagrange 余项。

欧拉方法

考虑一阶常微分方程的初值问题

前向欧拉方法

后向Euler算法或者后退Euler算法

梯形公式

改进的欧拉算法

其中，后向欧拉算法和梯形公式是隐式算法，前向欧拉算法和改进的欧拉算法是显式算法。欧拉算法计算容易，但是精度低，梯形公式精度高，但是是隐式形式，不易求解。将两式结合，则可以得到改进的欧拉公式。先用欧拉公式求出yn+1的一个粗糙的估计值，再用梯形方法进行精确化，称为校正值。

四阶龙格库塔算法（Runge-kutta method）

龙格-库塔算法是一种在工程上应用广泛的高精度单步算法，算法精度较高，如果预先取四个点就是四阶龙格-库塔算法。

欧拉方法和龙格-库塔算法的matlab实现

前向欧拉算法和后向欧拉算法

常微分方程y' = -y + x + 1;

x0 = 0;        %时间起始
x1 = 30;      %时间终止
h = 2.1;       %间隔
y0 = 1;    %初始值

代码：

%func.m文件  常微分方程为y' = -y + x + 1
function z = func(x, y)
z = -y + x + 1;

%forwardEuler.m文件
function [x, y] = forwardEuler(x0, y0, x1, h)
n = floor((x1 - x0)/h);
x = zeros(n + 1, 1);
y = zeros(n + 1, 1);
x(1) = x0;
y(1) = y0;
for i = 1 : n
    x(i + 1) = x(i) + h;
    y(i + 1) = y(i) + h * (-y(i) + x(i) + 1);
end

%backwardEuler.m文件    
function [x, y] = backwardEuler(x0, y0, x1, h)
n = floor((x1 - x0)/h);
x = zeros(n + 1, 1);
y = zeros(n + 1, 1);
x(1) = x0;
y(1) = y0;
for i = 1 : n
    x(i + 1) = x(i) + h;
    y(i + 1) = 1/(1 + h) * (y(i) + h * x(i + 1) + h);
end


%命令
x0 = 0;                                          %时间起始
x1 = 30;                                         %时间终止
h = 2.1;                                         %时间间隔
y0 = 1;                                          %y起始值
[x, y] = forwardEuler(x0, y0, x1, h);
[p, q] = backwardEuler(x0, y0, x1, h);
hold on                                          %将多个图均显示出来
plot(x, y,'r','LineWidth', 2)
plot(p, q, 'b','LineWidth', 2)
l = x0 : 0.01 : x1;
lu = exp(-l) + l;
plot(l, lu, 'g','LineWidth', 2)
legend('forwardEuler','backwardEuler','Theroy')

改进欧拉算法

常微分方程y' = y - 2x/y；

x0 = 0; 

x1 = 1；   

h = 0.1；

y0 = 1;

代码：

%常微分方程y' = y - 2x/y;

function z = func(x, y)
z = y - 2 * x / y;

useEuler.m文件  %向前欧拉算法代码
function [x, y] = useEuler(x0, x1,  y0, h)
n = floor((x1 - x0)/h);
x = zeros(n + 1, 1);
y = zeros(n + 1, 1);
x(1, 1) = x0;
y(1, 1) = y0;
for i = 2 : n+1
    x(i, 1) = x(i - 1, 1) + h;
    y(i, 1) = y(i - 1, 1) + h * (y(i - 1, 1) - 2 * x(i - 1, 1) / y(i - 1, 1));
end


improvedEuler.m  %改进的欧拉算法代码
function [x, y] = improvedEuler(x0, x1, y0, h)
n = floor((x1 - x0)/h);
x = zeros(n + 1, 1);
y = zeros(n + 1, 1);
x(1, 1) = x0;
y(1, 1) = y0;
for i = 2 : n+1
    x(i, 1) = x(i - 1, 1) + h;
    y(i, 1) = y(i - 1, 1) + h * (y(i - 1, 1) - 2 * x(i - 1, 1) / y(i - 1, 1));
    y(i, 1) = y(i - 1, 1) + h / 2 * (y(i - 1, 1) - 2 * x(i - 1, 1) / y(i - 1, 1) + y(i, 1) - 2 * x(i, 1) / y(i, 1));%用上一步求得的值代入这一步，变为显式算法
end

untitle.m文件
x0 = 0;
x1 = 1;
y0 = 1;
h = 0.1;
[x, y] = useEuler(x0, x1, y0, h);
yy = sqrt(2 * x + 1);    %精确解
[x, q] = improvedEuler(x0, x1, y0, h);
hold on
plot(x, y, 'b');
plot(x, yy, 'r');
plot(x, q, 'b.');

结果：

龙格-库塔算法

常微分方程y' = (x - y)/2;

x0 = 0;
x1 = 3;
y0 = 1;
h = 0.25;

代码：

function [x, y] = runge(x0, x1, y0, h)
n = (x1 - x0) / h;
x = zeros(n + 1);
y = zeros(n + 1);
x(1) = x0;
y(1) = y0;
for i = 1:n
    x(i + 1) = x(i) + h;
    k1 = fun(x(i), y(i));
    k2 = fun(x(i) + 0.5*h, y(i) + k1*h/2);
    k3 = fun(x(i) + 0.5*h, y(i) + k2*h/2);
    k4 = fun(x(i)+ h, y(i) + k3*h);
    y(i + 1) = y(i) + h*(k1 + 2*k2 + 2*k3 + k4)/6;
end
end


function z = fun(x, y)
z = (x - y)/2;
end


x0 = 0;
x1 = 3;
y0 = 1;
h = 0.25;
[x, y] = runge(x0, x1, y0, h);
[p, q] = ode45(@fun, [0, 3], y0);
hold on
plot(p, q, 'r');
plot(x, y, 'b.');

结果：

matlab中常用函数的总结

取整函数

floor函数朝负无穷方向取整，例如floor(-1.3) = -2, floor(2.8) = 2

ceil函数朝正无穷方向取整，例如ceil(1.2) = 2, ceil(-1.9) = -1

round函数 四舍五入, 例如round(3.3)  = 3, round(-4.8) = -5

fix函数朝0方向取整, 例如fix(1.8) = 1, fix(-1.8) = -1

绘图时常用的函数

plot函数plot函数的常用的调用语法如下:

plot(Y)，如果Y为实数向量，其维数为m，则plot(Y)等价于plot(X, Y), 其中X = 1:m。

例如，plot(t, y + 0.5)

plot(X1, Y1, X2, Y2.....)，Xi和Yi成对出现，该函数分别按照顺序取两个数据Xi和Yi绘图。

plot(X1, Y1, LineSpec,...)，将按顺序分别绘出由3个参数Xi、Yi和LineSpec定义的线条。其中参数LineSpec指明了线条的类型，再分别将配对的向量绘制出来。

例如，plot(x, y, axis('equal'), xlabel('x'), ylabel('y'))

曲线的线型和颜色参数

-, :, -., --,分别代表的是实线，虚线，点化线，双划线。

b, g, r, c, m, y, k, w,分别代表的是蓝，绿，红，青， 品红，黄，黑，白色。

hold on是当前轴及图形保持而不被刷新，准备接受此后的绘制。即在坐标轴内画了一个图后，再画另一个图时，原图还在，与新图共存；

hold off是使当前轴及图形不再具备被刷新的性质。即你在当前的坐标系中画了一幅图，此时的状态是hold off，则再画另一个图时，在轴上绘制的是新图，原图被替换了。

legend(S1, S2,...)绘制曲线所用线型、色彩或数据点型图例

subplot(m, n, p)或subplot(mnp), 子图函数所有的图按m行n列来排放，p表示该图所在的位置。