用Euler法解决初值问题

以如下的初值问题为例：
u′=4tu1/2,0≤t≤2u(0)=1\begin{aligned} &u^{'} = 4tu^{1/2},0\leq t \leq 2\\ &u(0)=1 \end{aligned} u′=4tu1/2,0≤t≤2u(0)=1
该方程的真实解为：
u(t)=(1+t2)2u(t)=(1+t^2)^2 u(t)=(1+t2)2
如下我们进行数值求解,将区间[0,2]进行分划，步长h=0.1则：
f(tn,un)=4nhun1/2f(t_{n},u_{n})=4nhu_n^{1/2} f(tn,un)=4nhun1/2

function [y] = f(u,t)
%u的导函数
y = 4*t*sqrt(u);
end

显式欧拉法

公式如下：
un+1=un+hf(tn,un)u_{n+1} =u_{n}+hf(t_{n},u_{n}) un+1=un+hf(tn,un)
代码如下：

function [x,y] = Euler(f,a,b,h,y0)
%Euler算法
%f为u的导函数，[a,b]为区间，h为步长，y0为初值N=(b-a)/h;x=a:h:b;y=zeros(1,N+1);y(1)=y0;for i=1:Ny(i+1)=y(i)+h*f(y(i),(i-1)*h);end
end

隐式欧拉法

公式如下：
un+1=un+hf(tn+1,un+1)u_{n+1} =u_{n}+hf(t_{n+1},u_{n+1}) un+1=un+hf(tn+1,un+1)
我们采用预估-校正的方法进行求解，以显格式作为预优算式( P )也即如下形式：
P:un+1[0]=un[M]+hfn[M]E:fn+1[m]=f(tn+1,un+k[m])C:un+1[m+1]=unM+hfn+1[M],0≤m≤M−1E:fn+1[M]=f(tn+1,un+k[M])\begin{aligned} &P:u_{n+1}^{[0]} =u_{n}^{[M]}+hf_{n}^{[M]}\\ &E:f_{n+1}^{[m]}=f(t_{n+1},u_{n+k}^{[m]})\\ &C:u_{n+1}^{[m+1]}=u_n^{M}+hf_{n+1}^{[M]},0\leq m \leq M-1\\ &E:f_{n+1}^{[M]}=f(t_{n+1},u_{n+k}^{[M]})\\ \end{aligned} P:un+1[0]=un[M]+hfn[M]E:fn+1[m]=f(tn+1,un+k[m])C:un+1[m+1]=unM+hfn+1[M],0≤m≤M−1E:fn+1[M]=f(tn+1,un+k[M])
代码如下：

function [x,y] = imp_Euler(f,a,b,h,y0,m)
%impEuler算法
% f为u的导函数，[a,b]为区间，h为步长，y0为初值，m为教证次数N=(b-a)/h;x=a:h:b;y=zeros(1,N+1);y(1)=y0;for i=1:N%Pif i==1y(i+1)=y(i)+h*f(y(i),(i-1)*h);elsey(i+1)=y(i)+h*E;endfor j=1:m%EE=f(y(i+1),i*h);%Cy(i+1)=y(i)+h*E;end%EE=f(y(i+1),i*h);end
end

改进Euler法（梯形公式）

公式如下：
un+1=un+h2[f(tn,un)+f(tn+1,un+1)]u_{n+1} =u_{n}+\frac{h}2[f(t_{n},u_{n})+f(t_{n+1},u_{n+1})] un+1=un+2h[f(tn,un)+f(tn+1,un+1)]

同样利用预估-校正算法求解
代码如下：

function [x,y] = improve(f,a,b,h,y0,m)
%改近Euler算法
%f为u的导函数，[a,b]为区间，h为步长，y0为初值，m为教证次数N=(b-a)/h;x=a:h:b;y=zeros(1,N+1);y(1)=y0;for i=1:N%Pif i==1y(i+1)=y(i)+0.5*h*f(y(i),(i-1)*h);elsey(i+1)=y(i)+0.5*h*E;endfor j=1:m%EE=f(y(i+1),i*h);%Cy(i+1)=y(i)+0.5*h*E+0.5*h*f(y(i),(i-1)*h);end%EE=f(y(i+1),i*h);end
end

main.m如下：

[x,y1]=Euler(@f,0,2,0.1,1);
[~,y2]=imp_Euler(@f,0,2,0.1,1,2);
[~,y3]=improve(@f,0,2,0.1,1,2);
g=@(z)(1+z.^2).^2;
z=x;
y=g(z);
figure(1)
plot(x,y1);
hold on
plot(x,y2)
plot(x,y3)
plot(x,y)
scatter(x,y1)
scatter(x,y2)
scatter(x,y3)
scatter(x,y,'*')
legend('E','impE','improve','gt','Ep','impEp','imporvep','gtp')
title('Euler法相关（隐式采用P（EC）^mE）')
hold off

结果图

参考文献

[1]李荣华，刘播. 微分方程数值解法（第四版）.北京：高等教育出版社，2009.

Euler法解微分方程相关推荐

龙格-库塔（Runge-Kutta）法解微分方程
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法.由于此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂.该算法是构建在数学支持的基础之上的. 对于一阶精度的欧 ...
Python06 向前Euler法、向后Euler法、梯形方法、改进的Euler方法以及四阶Runge_Kutta方法(附代码)
1. 实验结果 (1)解如下常微分方程: (2)分别使用向前 Euler 法.向后 Euler 法.梯形方法.改进的 Euler 方法以及四阶 Runge_Kutta 方法,结果如下图所示: 由结果 ...
龙格库塔(runge-kutta，RK)法求解微分方程
求解微分方程的意思就是,已知导函数,求原函数. 先声明一点,欧拉法.中值法.龙哥库塔法求解微分方程,得出的结果不是表达式,而是一系列离散点. 一.欧拉法递推问题描述:已知y'=f(x,y),求y(x ...
matlab 向前欧拉公式,向前欧拉公式在Matlab解微分方程初值解的问题
向前欧拉公式在Matlab解微分方程初值解的问题0 fuqilin1202013.07.04浏览527次分享举报用向前欧拉公式(10.8)求解初值问题,dy/dx=-3x+8x-7,y(0)=1,分 ...
2021-08-23 Euler法
数值分析--第一章:常微分方程初值问题数值方法第一章MATLAB与数值分析之常微分方程初值问题的数值方法. 本章内容主要参考的是李荣华老师的书 (参见[1]) 中内容第一章,感兴趣的读者可以参考该 ...
Matlab符号运算 - 解微分方程
什么是微分方程?如下: 如下图,也是一个微分方程: matlab使用 dsolve 命令解微分方程: 下面来解 dy/dx=4*x^3 这个方程: 一阶,二阶,多阶,都可以解:
matlab lu分解求线性方程组_计算方法（二）直接三角分解法解线性方程组
封面是WH2里春希在编辑部的上司麻理前辈,有一说一,这条线的第一次H有点恶趣味,不是很喜欢. 一:概述矩阵分解我学过的挺多种,比如极分解,谱分解,满秩分解,正交三角分解还有这里的直接三角分解大部分我 ...
matlab中函数或变量无法识别怎么办_用MATLAB巧解微分方程实例分析
点"考研竞赛数学"↑可每天"涨姿势"哦! MATLAB巧解微分方程实例分析王少华西安电子科技大学微分方程求解难, 字母一堆看着烦. 写错数字一时爽, 一直 ...
分枝定界法解0/1背包问题
分枝定界法解0/1背包问题关键词:分支定界.0-1背包分枝定界法简介分枝定界法按照树形结构来组织解空间,展开节点后,有两种策略: 策略一.把节点加入 FIFO 队列当中: 策略二.把节点加入到堆 ...
[Matlab科学计算] 四阶Runge-Kutta法解常微分方程
四阶Runge-Kutta法格式的详细推导请查找相关数值分析书籍,这里直接给出四阶Runge-Kutta法的经典格式和Matlab代码 Matlab代码如下:自行修改常微分方程即可 %% 四阶Rung ...

Euler法解微分方程

微分方程数值解法的一些练习