向前欧拉公式计算matlab程序函数
欧拉法(euler)求解常微分方程的matlab程序及案例
欧拉法(euler)求解常微分方程的matlab程序及案例欧拉方法是最初用于求解常微分方程的数值方法之一，它是一种显式的一步法，具有易于实施的优点，特别适合初学者使用。
本文将介绍欧拉法的原理和使用MATLAB求解常微分方程的具体方法，同时给出一个简单的实例进行说明。
一、欧拉法原理考虑一个一阶常微分方程y'=f(t,y)，欧拉法的基本思想是将时间步长Δt均分成n个小步长，从y(t0)开始依次计算每个时刻的值，得到一列估计值y1, y2, …, yn。
欧拉法的计算公式为：（1）y1=y(t0+Δt)=y(t0)+Δtf(t0, y0)（2）y2=y(t0+2Δt)=y(t0+Δt)+Δtf(t0+Δt, y1)（3）yn=y(t0+nΔt)=y(t0+(n-1)Δt)+Δtf(t0+(n-1)Δt, yn-1)可以看出，欧拉法的核心在于利用已知的t和y计算f(t,y)，从而获得y的逼近值。
但是需要注意的是，步长Δt越小，计算所需的时间和内存就越多，而精度却并不一定提高。
因此在实际应用中需要结合具体问题选择合适的步长。
二、MATLAB求解常微分方程的具体方法（1）定义常微分方程我们以一个简单的例子开始，考虑求解y'=1-y，y(0)=0.5在[0,1]区间内的积分。
首先定义匿名函数dydt，将其传到ode45中求解：dydt=@(t,y)1-y;[t,y]=ode45(dydt,[0 1],0.5);plot(t,y,'-o')运行以上代码可以得到结果，其中plot函数用于绘制图像。
但是，由于求解过程中计算机执行到ode45函数时可能需要很长时间，因此需要更快捷的方法。
（2）利用欧拉法求解方程欧拉法求解方程首先需要定义步长Δt，这里设Δt为0.1。
定义起始值y=[0.5]、时间向量t=0:Δt:1，然后计算列向量y的估计值：t=0:0.1:1;y=zeros(size(t));y(1)=0.5;for n=1:length(t)-1y(n+1)=y(n)+0.1*(1-y(n));endplot(t,y,'-o')以上代码的执行结果与前面的ode45方法相同，但是速度更快。
用MATLAB程序生动地演示欧拉公式
⽤MATLAB程序⽣动地演⽰欧拉公式下⾯的MA TLAB 程序⽣动地演⽰欧拉公式Exp(t) = cos(t) + j sin(t)% Henry-104% 本程序演⽰欧拉公式% Jan.25th,2012%h_fig1 = figure;set(h_fig1, 'unit', 'normalized', 'position', [0.1, 0.1, 0.9, 0.9]);set(h_fig1, 'defaultuicontrolunits', 'normalized');h_text1 = uicontrol(h_fig1, 'Style', 'text', 'Position', [0.71, 0.73, 0.25, 0.05],... % 创建⽂本框'String', '▲是cos 曲线的起点', 'ForegroundColor', 'r', 'FontName', '⿊体',...'FontSize', 12, 'FontWeight', 'Bold', 'BackgroundColor', [1, 1, 1]);h_text2 = uicontrol(h_fig1, 'Style', 'text', 'Position', [0.71, 0.78, 0.25, 0.05],... % 创建⽂本框'String', 'Δ是sin 和exp 曲线的起点', 'ForegroundColor', 'r', 'FontName', '⿊体',...'FontSize', 12, 'FontWeight', 'Bold', 'BackgroundColor', [1, 1, 1]);h_pushbutton1 = uicontrol(h_fig1, 'Style', 'PushButton', 'Position', [0.82, 0.12, 0.07, 0.06],...'string', '退出', 'BackgroundColor', [0.8 0.9 0.8], 'ForegroundColor', 'r', 'FontSize', 14, 'FontWeight', 'Bold',...'callback', 'delete(h_fig1),')h_axes0 = axes('Box', 'on', 'Position', [0.15, 0.18, 0.56, 0.68], 'FontSize', 8)set(gcf,'color','w');w = 0.1*pit = 0:40; % 在前进⽅向绕了2 圈%a = -ones(1,length(t));plot3(cos(w*t),t,sin(w*t),'b', 'LineWidth', 2);grid on; hold on;hc = plot3(cos(w*t),t,a,'k--'); hold on;set(hc, 'Color', 'r', 'LineWidth', 2);a=-a;hs = plot3(a,t,sin(w*t),'r-.'); hold on;set(hs, 'Color', 'k', 'LineWidth', 2);text(0.7,0.3,0.6, ' <-- CCW', 'FontSize', 14, 'FontWeight', 'Bold'); text(1,0,-1, ' ▲Cos', 'Color', 'r', 'FontSize', 14, 'FontWeight', 'Bold'); text(1,0,0, ' Δ Sin', 'FontSize', 14, 'FontWeight', 'Bold');%xlabel('x', 'FontSize', 14, 'FontWeight', 'Bold');ylabel('t', 'FontSize', 14, 'FontWeight', 'Bold');zlabel('y', 'FontSize', 14, 'FontWeight', 'Bold');title('演⽰欧拉公式y = exp(jwt) = cos(wt) + jsin(wt)', 'Color', 'b', …'FontSize', 18, 'FontWeight', 'Bold');%line([-1,-1],[39.9,39.9],[-1,1],'LineWidth',3, 'Color', 'r');line([1,1],[39.9,39.9],[-1,1],'LineWidth',3, 'Color', 'r');line([-1,-1],[0,0],[-1,1],'LineWidth',3, 'Color', 'r');line([1,1],[0,0],[-1,1],'LineWidth',3, 'Color', 'r');line([-1,-1],[0,40],[-1,-1],'LineWidth',3, 'Color', 'k');line([-1,1],[0,0],[-1,-1],'LineWidth',3, 'Color', 'b')line([-1,1],[40,40],[1,1],'LineWidth',3, 'Color', 'b')line([-1,1],[40,40],[-1,-1],'LineWidth',3, 'Color', 'b')line([-1,1],[0,0],[1,1],'LineWidth',3, 'Color', 'b')line([-1,1],[0,0],[0,0],'LineWidth',2, 'Color', 'k');line([0,0],[0,0],[-1,1],'LineWidth',2, 'Color', 'k');line([0,0],[40,40],[-1,1],'LineWidth',2, 'Color', 'k');line([0,0],[0,40],[0,0],'LineWidth',2, 'Color', 'k');line([-1,1],[40,40],[0,0],'LineWidth',2, 'Color', 'k');line([0,0],[0,40],[0,0],'LineWidth',2, 'Color', 'k');text(0,0,0.12,'O', 'FontSize', 14, 'FontWeight', 'Bold', 'Color', 'r') text(0,40,0.12,'O', 'FontSize', 14, 'FontWeight', 'Bold', 'Color', 'b')程序运⾏结果如下所⽰。
欧拉法求解一阶微分方程matlab
为了更好地理解欧拉法求解一阶微分方程在Matlab中的应用，我们首先来了解一些背景知识。
一阶微分方程是指只含有一阶导数的方程，通常表示为dy/dx=f(x,y)，其中f(x,y)是关于x和y的函数。
欧拉法是一种常见的数值解法，用于求解微分方程的近似数值解。
它是一种基本的显式数值积分方法，通过将微分方程转化为差分方程来进行逼近。
在Matlab中，我们可以利用欧拉法求解一阶微分方程。
我们需要定义微分方程的函数表达式，然后选择合适的步长和初始条件，最后使用循环计算逼近解。
下面我们来具体讨论如何在Matlab中使用欧拉法来求解一阶微分方程。
我们假设要求解的微分方程为dy/dx=-2x+y，初始条件为y(0)=1。
我们可以通过以下步骤来实现：1. 我们需要在Matlab中定义微分方程的函数表达式。
在Matlab中，我们可以使用function关键字来定义函数。
在这个例子中，我们可以定义一个名为diff_eqn的函数，表示微分方程的右侧表达式。
在Matlab中，这个函数可以定义为：```matlabfunction dydx = diff_eqn(x, y)dydx = -2*x + y;end```2. 我们需要选择合适的步长和初始条件。
在欧拉法中，步长的选择对于数值解的精度非常重要。
通常情况下，可以先尝试较小的步长，然后根据需要进行调整。
在这个例子中，我们可以选择步长h=0.1，并设置初始条件x0=0，y0=1。
3. 接下来，我们可以使用循环来逼近微分方程的数值解。
在每一步，根据欧拉法的迭代公式y(i+1) = y(i) + h * f(x(i), y(i))，我们可以按照下面的Matlab代码计算逼近解：```matlabh = 0.1; % 步长x = 0:h:2; % 定义计算区间y = zeros(1, length(x)); % 初始化y的值y(1) = 1; % 设置初始条件for i = 1:(length(x)-1) % 欧拉法迭代y(i+1) = y(i) + h * diff_eqn(x(i), y(i));end```通过上述步骤，在Matlab中就可以用欧拉法求解一阶微分方程。
欧拉方法及其改进的欧拉方法的Matlab实现
11( n n y x y ++−。为了估计它,由Taylor展开得到的精确值1( n y x +是
2'
''
31( ( ( ( ( 2
n n n n h y x y x hy x y x O h +=+++ (5
2.欧拉方法、改进的欧拉方法及Matlab实现
下面主要讨论一阶常微分方程的初值问题,其一般形式为:
' 00
(,
( y f x y y x y ⎧=⎨
=⎩ (1我们知道,只要函数(, f x y适当光滑——譬如关于y满足利普希茨(Lipschitz条件
(, (, f x y f x y L y y −≤−
改进的欧拉方法是先用欧拉公式求1( n y x +的一个近似值1n y +,称为预测值,然后用梯形公式进行矫正并求得近似值1n y +。即
1111(, [(, (, ]
2n n n n n n n n n n y y f x y h h
y y f x y f x y ++++⎧=+⎪
⎨=++⎪⎩
(8 2.2.2改进的欧拉方法的误差估计
方法是一阶方法,因此它的精度不高。
2.2改进的欧拉方法
2.2.1改进的欧拉方法
用数值积分方法离散化问题(1,两端积分可得
1
1( ( (, ( (0,1, 2, n n
x n n x y x y x f x y x dx n ++−==∫
欧拉法
欧拉法的MATLAB实现班级：123082--04 姓名：应业敏学号：20081001088摘要：在数学和计算机科学中，欧拉方法（Euler method）命名自它的发明者莱昂哈德·欧拉，是一种一阶数值方法，用以对给定初值的常微分方程(即初值问题)求解。
它是一种解决常微分方程数值积分的最基本的一类显型方法(Explicit method)。
欧拉法（euler method）是以流体质点流经流场中各空间点的运动即以流场作为描述对象研究流动的方法。
它不直接追究质点的运动过程，而是以充满运动液体质点的空间——流场为对象。
研究各时刻质点在流场中的变化规律。
将个别流体质点运动过程置之不理，而固守于流场各空间点。
通过观察在流动空间中的每一个空间点上运动要素随时间的变化，把足够多的空间点综合起来而得出的整个流体的运动情况。
基本思想：基本思想是迭代。
其中分为前进的EULER法、后退的EULER法、改进的EULER法。
所谓迭代，就是逐次替代，最后求出所要求的解，并达到一定的精度。
误差可以很容易的计算出来。
特点：单步，显式，一阶求导精度，截断误差贰阶缺点：欧拉法简单地取切线的端点作为下一步的起点进行计算，当步数增多时，误差会因积累而越来越大。
因此欧拉格式一般不用于实际计算。
改进欧拉格式：为提高精度，需要在欧拉格式的基础上进行改进。
采用区间两端的函数值的平均值作为直线方程的斜率。
改进欧拉法的精度为二阶。
算法为：微分方程的本质特征是方程中含有导数项，数值解法的第一步就是设法消除其导数值，这个过程称为离散化。
实现离散化的基本途径是用向前差商来近似代替导数，这就是欧拉算法实现的依据。
欧拉(Euler)算法是数值求解中最基本、最简单的方法，但其求解精度较低，一般不在工程中单独进行运算。
所谓数值求解，就是求问题的解y(x)在一系列点上的值y(xi)的近似值yi。
对于常微分方程：，x∈[a,b]y(a) = y0可以将区间[a,b]分成n段，那么方程在第x i点有y'(x i) = f(x i,y(x i))，再用向前差商近似代替导数则为：，在这里，h是步长，即相邻两个结点间的距离。
matlab实例讲解欧拉法求解微分方程
欧拉法是数值分析中常用的一种方法，用于求解常微分方程的数值解。
在MATLAB中，可以通过编写相应的代码来实现欧拉法求解微分方程。
下面我们将通过具体的实例来讲解MATLAB中如何使用欧拉法求解微分方程。
我们要了解欧拉法的基本原理。
欧拉法是一种通过迭代逼近微分方程解的方法，它基于微分方程的定义，通过离散化的方法逼近微分方程的解。
其基本思想是利用微分方程的导数定义，将微分方程以差分形式进行逼近。
具体而言，欧拉法通过将微分方程转化为差分方程的形式，然后通过迭代逼近得到微分方程的数值解。
接下来，我们通过一个具体的实例来讲解MATLAB中如何使用欧拉法求解微分方程。
假设我们要求解以下的一阶常微分方程：(1) dy/dx = x + y(2) y(0) = 1现在我们来编写MATLAB代码来实现欧拉法求解这个微分方程。
我们需要确定微分方程的迭代步长和迭代范围。
假设我们将x的范围取为0到10，步长为0.1。
接下来，我们可以编写MATLAB代码如下：```matlab欧拉法求解微分方程 dy/dx = x + y定义迭代步长和范围h = 0.1;x = 0:h:10;初始化y值y = zeros(1,length(x));y(1) = 1;使用欧拉法迭代求解for i = 1:(length(x)-1)y(i+1) = y(i) + h * (x(i) + y(i));end绘制图像plot(x,y,'-o');xlabel('x');ylabel('y');title('欧拉法求解微分方程 dy/dx = x + y');```在这段MATLAB代码中，我们首先定义了迭代的步长和范围，并初始化了微分方程的初始值y(0) = 1。
然后通过for循环使用欧拉法进行迭代求解微分方程，最后绘制出了微分方程的数值解的图像。
通过以上的实例讲解，我们可以看到，在MATLAB中使用欧拉法求解微分方程是非常简单而直观的。
matlab软件欧拉算法教程
y( xn1 ) y( xn ) hK
*
寻求计算平均斜率的算法
考察欧拉法，以xn的斜率值
K1 f ( xn , yn )
作为平均斜率
考察改进的欧拉法，可以将其改写为：
1 1 yi 1 yn h K1 K 2 2 2 K1 f ( x n , yn ) K 2 f ( xn h, yn hK1 )
Step 1: 将 K2 在 ( xn , yn ) 点作 Taylor 展开
K 2 f ( xn ph , yn phK1 ) f ( xn , yn ) phf x ( xn , yn ) phK1 f y ( xn , yn ) O( h2 )
y( xn ) phy( xn ) O(h2 )
Euler法
xn Yn
0.1 1.1 0.2 1.1918
改进Euler法
0.0046 1.0959 0.0086 1.1841 0.0005 0.0009
准确解
1.0954 1.1832
|yn-y(xn)| Yn
|yn-y(xn)| y(xn)
0.3 1.2774
0.4 1.3582 0.5 1.4351 0.6 1.5090 0.7 1.5803 0.8 1.6498 0.9 1.7178 1.0 1.7848
Y
x1 b
N 结束
例题3
例2
求解初值问题 步长h 0.1) ( 2x y y y y ( 0) 1 (0 x 1)
解
f ( x, y) y 2 x / y
y p yn hf ( xn , yn ) 改进Euler格式 yc yn hf ( xn1 , y p ) yn1 1 ( y p yc ) 2
matlab欧拉法求解微分方程
matlab欧拉法求解微分方程欧拉法是一种用来求解微分方程数值解的方法，它是由欧拉在18世纪提出的。
该方法基于微分方程的定义，将微分方程转化为差分方程，从而通过求解差分方程获得微分方程的数值解。
本文将介绍欧拉法的基本原理和实现步骤，并通过一个具体的例子来演示其应用。
欧拉法的基本原理是将微分方程中的导数近似为差商，从而将微分方程转化为差分方程。
对于一阶微分方程y'(x)=f(x,y(x))，我们可以将其转化为差分方程如下：y(x+h)≈y(x)+h*f(x,y(x)),其中h为步长，x为自变量，y为因变量，f为给定的函数。
根据该差分方程，我们可以通过递归的方式求解微分方程的数值解。
具体的求解步骤如下：1.确定微分方程的初始条件，即给定初始点(x0,y0)。
2.选择一个适当的步长h。
3.通过欧拉法的递推关系式计算y的近似值：y(n+1)=y(n)+h*f(x(n),y(n)),其中n表示当前迭代的步数，x(n)和y(n)分别表示第n步迭代的自变量和因变量的数值。
4.重复步骤3，直到达到所需的计算精度或结束条件。
下面，我们通过一个具体的例子来演示欧拉法的应用。
假设我们需要求解微分方程y'(x)=y(x)-x的数值解，初始条件为y(0)=1、我们可以利用欧拉法来计算该微分方程在区间[0,1]上的数值解。
首先，我们定义微分方程的函数f(x,y)=y-x。
然后，选择一个适当的步长h，比如h=0.1、根据欧拉法的递推关系式，我们可以得到下面的迭代公式：y(n+1)=y(n)+h*(y(n)-x(n)),其中n表示当前迭代的步数，x(n)和y(n)分别表示第n步迭代的自变量和因变量的数值。
接下来，我们可以利用该迭代公式来求解微分方程的数值解。
从初始点(x0,y0)=(0,1)开始，进行迭代计算。
具体的迭代过程如下：步骤1：设置初始点(x0,y0)=(0,1)。
步骤2：选择步长h=0.1步骤3：利用迭代公式计算数值解y(n+1)=y(n)+h*(y(n)-x(n))。
欧拉法(euler)求解常微分方程的matlab程序及案例
欧拉法(euler)求解常微分方程的matlab程序及案例欧拉法是一种常见的求解常微分方程的数值解法，在MATLAB中可以通过编写简单的程序实现。
本文将介绍欧拉法的MATLAB程序及应用案例。
首先，让我们考虑以下的常微分方程：dy/dx = f(x, y)其中y是关于x的函数，f是已知的函数。
我们可以通过欧拉法求解该方程。
欧拉法的基本思想是将区间[x0, xn]分成n等份，然后用以下式子计算y的值：y(i+1) = y(i) + h*f(x(i), y(i))其中h是步长，x(i)和y(i)分别表示当前的x和y值，y(i+1)表示下一个y值。
通过重复上述计算，欧拉法可以求出y在x=n处的值。
下面是欧拉法的MATLAB程序：% 默认参数x0 = 0; % 初始值xn = 1; % 终止值y0 = 1; % 初始y值h = 0.1; % 步长f = @(x, y) -y; % 函数n = (xn - x0) / h; % 时间步数x = x0; % 初始x值y = y0; % 初始y值for i = 1:ny = y + h * f(x, y);x = x + h;enddisp(['y在x = ', num2str(xn), '处的值为：',num2str(y)]);在上述程序中，我们定义了默认的初始值、终止值、初始y值和函数。
程序中的n表示时间步数，x和y分别表示当前的x和y值。
通过for循环，欧拉法可以重复计算y的值，并最终求出y在x=n处的值。
下面是一个用欧拉法求解dy/dx = -y的应用案例：% 默认参数x0 = 0; % 初始值xn = 5; % 终止值y0 = 1; % 初始y值h = 0.1; % 步长f = @(x, y) -y; % 函数n = (xn - x0) / h; % 时间步数x = x0; % 初始x值y = y0; % 初始y值% 初始化结果数组result = zeros(n + 1, 2);result(1,:) = [x0 y0];for i = 1:ny = y + h * f(x, y);x = x + h;% 保存结果result(i + 1,:) = [x y];end% 绘制图形plot(result(:,1), result(:,2), '-o');xlabel('x');ylabel('y');title('欧拉法求解dy/dx=-y');在上述案例中，我们使用默认的参数，求解dy/dx=-y的方程。
向前欧拉公式计算matlab程序函数
向前欧拉公式计算matlab程序函数向前欧拉公式计算matlab程序函数function [H,X,Y,k,h,P]=QEuler(funfcn,x0,b,y0,tol) %初始化.pow=1/3;if nargin < 5 | isempty(tol),tol = 1.e-6;end;x=x0;h=0.0078125*(b-x);y=y0(:);p=128;H=zeros(p,1);X=zeros(p,1);Y=zeros(p,length(y));k=1;X(k)=x;Y(k,:)=y';% 绘图.clc,x,h,y% end% 主循环.while (x<b)&(x+h>x)</b)&(x+h>if x+h>bh=b-x;end%计算斜率.fxy=feval(funfcn,x,y);fxy=fxy(:);%计算误差，设定可接受误差.delta=norm(h*fxy,'inf');wucha=tol*max(norm(y,'inf'),1.0);% 当误差可接受时重写解.if delta<=wuchax=x+h; y=y+h*fxy; k=k+1;if k>length(X)X=[X;zeros(p,1)];Y=[Y;zeros(p,length(y))];H=[H;zeros(p,1)];endH(k)=h;X(k)=x;Y(k,:)=y';plot(X,Y,'rp'),gridxlabel('自变量X'),ylabel('因变量Y')end%更新步长.if delta~=0.0h=min(h*8,0.9*h*(wucha/delta)^pow); endendif (x<b)< p="">disp('Singularity likely.'), xendH=H(1:k);X=X(1:k);Y=Y(1:k,:);n=1:k;P=[n',H,X,Y]endfunction y=funfcn(x,y)y(1)=-3*y+8*x-7;endclc,clear ,close allsubplot(2,1,1)x0=0;y0=1;b=2;n=10;h=2/10;[k,X,Y,P,REn]=QEuler(@funfcn,x0,y0,b,h)hold onS1=1+1/6*(6+12*X+30*exp(2*X)).^(1/2)plot(X,S1,'b-')title('用向前欧拉公式计算dy/dx=y-x/(3y)，y(0)=1在[0,2]上的数值解') legend('n=10数值解','精确解')hold offjdwucY=S1-Y;jwY=S1-Y;xwY=jwY./Y;k1=1:n;k=[0,k1];% P1=[k',X,Y,S1,jwY,xwY]subplot(2,1,2)n1=100;h1=2/100;[k,X1,Y1,P1,Ren1]=QEuler(@funfcn,x0,y0,b,h1) hold onS2=1+1/6*(6+12*X1+30*exp(2*X1)).^(1/2) plot(X1,S2,'b-') legend('n=100数值解','精确解')hold offjwY1=S2-Y1;xwY1=jwY1./Y1;k1=1:n1;k=[0,k1];% P2=[k',X1,Y1,S2,jwY1,xwY1]</b)<>。
欧拉函数及代码实现
欧拉函数及代码实现对正整数n，欧拉函数是⼩于n的正整数中与n互质的数的数⽬，它⼜称为Euler’s totient function、φ函数、欧拉商数等。
例如φ(8)=4，因为1,3,5,7均和8互质，特殊的（φ(1)=1）。
1.根据欧拉函数公式：euler(x) = x*(1-1/p1)(1-1/p2)……(1-1/pn),p为x的质因数。
⽤公式法求解代码int euler(int n){int res=n,a=n;for(int i=2;i*i<=a;i++){if(a%i==0){res=res/i*(i-1);//先进⾏除法是为了防⽌中间数据的溢出while(a%i==0) a/=i;}}if(a>1) res=res/a*(a-1);return res;}2.O(n^2)的筛法：#include<cstdio>const int N = 100000 + 5;int phi[N];void Euler(){phi[1] = 1;for(int i = 2; i < N; i ++){if(!phi[i]){for(int j = i; j < N; j += i){if(!phi[j]) phi[j] = j;phi[j] = phi[j] / i * (i-1);}}}}int main(){Euler();}3.线性筛法（最快）：#include<iostream>#include<cstdio>#define N 40000using namespace std;int n;int phi[N+10],prime[N+10],tot,ans;bool mark[N+10];void getphi(){int i,j;phi[1]=1;for(i=2;i<=N;i++)//相当于分解质因式的逆过程{if(!mark[i]){prime[++tot]=i;//筛素数的时候⾸先会判断i是否是素数。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
向前欧拉公式计算matlab程序函数
function [H,X,Y,k,h,P]=QEuler(funfcn,x0,b,y0,tol)
%初始化.
pow=1/3;
if nargin < 5 | isempty(tol),
tol = 1.e-6;
end;
x=x0;
h=0.0078125*(b-x);
y=y0(:);
p=128;
H=zeros(p,1);
X=zeros(p,1);
Y=zeros(p,length(y));
k=1;
X(k)=x;
Y(k,:)=y';
% 绘图.
clc,x,h,y
% end
% 主循环.
while (x<b)&(x+h>x)
if x+h>b
h=b-x;
end
%计算斜率.
fxy=feval(funfcn,x,y);
fxy=fxy(:);
%计算误差，设定可接受误差.
delta=norm(h*fxy,'inf');
wucha=tol*max(norm(y,'inf'),1.0);
% 当误差可接受时重写解.
if delta<=wucha
x=x+h; y=y+h*fxy; k=k+1;
if k>length(X)
X=[X;zeros(p,1)];
Y=[Y;zeros(p,length(y))];
H=[H;zeros(p,1)];
end
H(k)=h;
X(k)=x;
Y(k,:)=y';
plot(X,Y,'rp'),
grid
xlabel('自变量X'),
ylabel('因变量Y')
end
%更新步长.
if delta~=0.0
h=min(h*8,0.9*h*(wucha/delta)^pow);
end
end
if (x<b)
disp('Singularity likely.'), x
end
H=H(1:k);
X=X(1:k);
Y=Y(1:k,:);
n=1:k;
P=[n',H,X,Y]
end
function y=funfcn(x,y)
y(1)=-3*y+8*x-7;
end
clc,clear ,close all
subplot(2,1,1)
x0=0;
y0=1;
b=2;
n=10;
h=2/10;
[k,X,Y,P,REn]=QEuler(@funfcn,x0,y0,b,h)
hold on
S1=1+1/6*(6+12*X+30*exp(2*X)).^(1/2)
plot(X,S1,'b-')
title('用向前欧拉公式计算dy/dx=y-x/(3y)，y(0)=1在[0,2]上的数值解') legend('n=10数值解','精确解')
hold off
jdwucY=S1-Y;
jwY=S1-Y;
xwY=jwY./Y;
k1=1:n;
k=[0,k1];
% P1=[k',X,Y,S1,jwY,xwY]
subplot(2,1,2)
n1=100;
h1=2/100;
[k,X1,Y1,P1,Ren1]=QEuler(@funfcn,x0,y0,b,h1) hold on
S2=1+1/6*(6+12*X1+30*exp(2*X1)).^(1/2) plot(X1,S2,'b-')
legend('n=100数值解','精确解')
hold off
jwY1=S2-Y1;
xwY1=jwY1./Y1;
k1=1:n1;k=[0,k1];
% P2=[k',X1,Y1,S2,jwY1,xwY1]。