最近参加了数学建模,对于老师说的Euler算法的不同步长的精度不一样,编写了一个M函数文件来实现这个精度的比较,把函数附上:
function [x,y]= Euler(varargin)
%这里使用可变输出输入函数的
%varargin{1}为求解常微分方程的表达式
%varargin{2}为求解常微分方程的定解条件
%需要给出的变量有常微分方程的范围a,b(varargin{3},varargin{4})
%n为对这个区间的分割(varargin{5})
%xlt写于7月8日
%取得算法需要的变量,并附上容易理解的含义变量
a = varargin{3};
b = varargin{4};%自变量的范围
n = varargin{5};%区间的分割次数
h = (b - a)/n;%步长
Dy = varargin{1}; %常微分方程的表达式
y0 = varargin{2}; %常微分方程的定解条件表达式
%首先求出所给常微分方程问题的精确解
x1 = zeros(n+1,1);
y1 = zeros(n+1,1);
syms f1; syms x;
f1 = dsolve(Dy,y0,'x');
x1(1) = a;
y1(1) = subs(f1,{x},{x1(1)});
for i = 2:(n+1)
x1(i) =
x1(i-1) + h;
y1(i) =
double(subs(f1,{x},{x1(i)}));
end
%利用Euler方法求解近似数值微分解
x2 = zeros(n+1,1);
y2 = zeros(n+1,1);
syms y;
x2(1) = a;
y2(1) = subs(f1,{x},{a});%获得原方程的初解
for i = 2:(n+1)
x2(i) =
x2(i-1) + h;
y2(i) =
y2(i-1) + h .*
double(subs(Dy(5:end),{x,y},{x2(i-1),y2(i-1)}));%特别记录Matlab中的字符串操作,提取子字符串即A(3:6)...
end
%返回经过Euler算法算出x与y的值
x = x2;
y = y2;
%画图进行误差比较
plot(x1,y1,'r');
hold on;
plot(x2,y2,'b');
特此记录,以后写了新的算法再分享