| matlab分别用欧拉法，改进欧拉法，和四阶龙格法对同一个函数进行拟合 | 您所在的位置：网站首页 › 改进euler公式 › matlab分别用欧拉法，改进欧拉法，和四阶龙格法对同一个函数进行拟合 |
matlab分别用欧拉法，改进欧拉法，和四阶龙格法对同一个函数进行拟合
|
由于数值分析临时布置了有关于MATLAB一个大作业，才发现身边有相当一部分的同学对matlab相当不熟悉，于是做了一个关于这三个方法的小文章做一点抛砖引玉的作用。 首先是关于这三个方法的总结，欧拉法属于显式的算法，即可以直接由已知的数据推出所需要的数据，但是贴合的程度较低。改进欧拉法增加了一步计算量，使曲线能够更加的贴近于原数据，而四阶龙格法则是贴合程度最高，但相对计算量最大的一个算法。 各个算法的公式如下： 一阶微分方程：
欧拉法：
改进欧拉法：
四阶龙格法：
下面直接上MATLAB代码来实现操作，具体操作的步骤和原因看代码的注释， %func1.m文件 常微分方程为y' = 2*y/(x+1) function z = func1(x, y) z = 2*y/(x+1); %欧拉法函数定义 %forwardEuler1.m文件 向前欧拉法 function [x, y] = forwardEuler(x0, y0, x1, h) %定义函数 n = floor((x1 - x0)/h); %得到所要建立的向量长度 x = zeros(n + 1, 1); %建立x和y的向量 y = zeros(n + 1, 1); x(1) = x0; %赋值：由于matlab是从1开始而不是从0开始 y(1) = y0; % 带入公式，用for循环得出向量各个值 for i = 1 : n x(i + 1) = x(i) + h; y(i + 1) = y(i) + h * (2 * y(i) / (x(i)+1)); end %命令1 x0 = 0; %x起始值 x1 = 1; %x终止值 h = 0.1; %步长 y0 = 1; %y起始值 [x, y] = forwardEuler(x0, y0, x1, h); %函数代入 hold on %将多个图均显示出来 plot(x, y,'r','LineWidth', 2) %画出图像，红色，线的宽度为2 l = x0 : 0.1 : x1; %取标准值 lu = (x+1).*(x+1); plot(l, lu, 'g','LineWidth', 2) legend('欧拉方法','理论值') %func2.m文件 常微分方程为y' = 2*y/(x+1) function z = func2(x, y) z = 2*y/(x+1); %改进欧拉法函数定义 %Euler2.m文件 未改进的欧拉法 function [x, y] = Euler(x0, x1, y0, h) n = floor((x1 - x0)/h); % 求n值 x = zeros(n + 1, 1); % 定义x与y向量 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 * (2 * y(i - 1, 1) / (x(i - 1, 1) + 1)); end %improvedEuler2.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) = x0; y(1) = y0; % 改进欧拉法公式代入，由于有不同的公式则先得出z1，才能进行z2的运算 for i = 1 : n x(i + 1) = x(i) + h; z1 = 2 * y(i) / (x(i)+1); z2 = 2 * (y(i) + h * z1) / (x(i) + h + 1); y(i + 1) = y(i) + h / 2 * (z1 + z2); end %命令2 x0 = 0; x1 = 1; y0 = 1; h = 0.1; [x, y] = useEuler(x0, x1, y0, h); yy = (1+x).*(1+x); %精确解 [x, q] = improvedEuler(x0, x1, y0, h); hold on plot(x, y, 'k'); plot(x, yy, 'r'); plot(x, q, 'b'); legend('传统欧拉法','理论值','改进欧拉法') %func3.m文件 常微分方程为y' = 2*y/(x+1) function z = func3(x, y) z = 2*y/(x+1); %四阶龙格法函数定义 %runge3.m文件 四阶龙格法方法展示 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 = func(x(i), y(i)); k2 = func(x(i) + 0.5*h, y(i) + k1*h/2); k3 = func(x(i) + 0.5*h, y(i) + k2*h/2); k4 = func(x(i)+ h, y(i) + k3*h); y(i + 1) = y(i) + h*(k1 + 2*k2 + 2*k3 + k4)/6; end end %命令3 x0 = 0; x1 = 1; y0 = 1; h = 0.1; [x, y] = runge(x0, x1, y0, h); yy = (1+x).*(1+x); %精确解 hold on plot(x, yy, 'r'); plot(x, y, 'b.'); %由于四阶龙格法与理论值的贴合度太高，所以一个采用了线段形式，一个采用了点的形式 %其中点为四阶龙格，线是标准值其中根据1.2.3分别进行不同的分文件选取，即后缀含有1的代表欧拉法展示，后缀含有2的代表改进欧拉法的展示，后缀带3的代表四阶龙格法展示，大概会产生十个分文件，分别在命令的文件中进行运行即可得到所需的下列三个图像。
|
【本文地址】
公司简介
联系我们