Matlab求解边值问题方法 例题解读
重要:MATLAB常微分方程(组)数值解法
![重要:MATLAB常微分方程(组)数值解法](https://img.taocdn.com/s3/m/713c0de9284ac850ac024246.png)
Matlab常微分方程求解问题分类
边值问题:
初值问题:
• 定解附加条件在自变量 的一端
• 一般形式为: y' f (x, y)
y(a)
y0
• 初值问题的数值解法一 般采用步进法,如 Runge-Kutta法
➢ 在自变量两端均给定附加 条件
y' f (x, y)
➢ 一般形式:y(a)y1, y(b)y2
1.根据常微分方程要求的求解精度与速度要求
求解初值问题:
y
'
y
2x y
y ( 0 ) 1
(0x1)
比较ode45和ode23的求解精度和速度
ode45和ode23的比较-1
function xODE clear all clc
format long
y0 = 1; [x1,y1] = ode45(@f,[0,1],y0); [x2,y2] = ode23(@f,[0,1],y0); plot(x1,y1,'k-',x2,y2,'b--') xlabel('x') ylabel('y')
rD = k(3)*C(2)-k(5)*C(4);
rE = k(4)*C(3)+k(5)*C(4);
% Mass balances dCdt = [rA; rB; rC; rD; rE];
三个串联的CSTR等温反应器(例4-3)
function IsothermCSTRs clear all clc CA0 = 1.8; % kmol/m^3 CA10 = 0.4; % kmol/m^3 CA20 = 0.2; % kmol/m^3 CA30 = 0.1; % kmol/m^3 k = 0.5; % 1/min tau = 2; stoptime = 2.9; % min [t,y] = ode45(@Equations,[0 stoptime],[CA10 CA20 CA30],[],k,CA0,tau); disp(' Results:') disp(' t CA1 CA2 CA3') disp([t,y]) plot(t,y(:,1),'k--',t,y(:,2),'b:',t,y(:,3),'r-') legend('CA_1','CA_2','CA_3') xlabel('Time (min)') ylabel('Concentration') % -----------------------------------------------------------------function dydt = Equations(t,y,k,CA0,tau) CA1 = y(1); CA2 = y(2); CA3 = y(3); dCA1dt = (CA0-CA1)/tau - k*CA1; dCA2dt = (CA1-CA2)/tau - k*CA2; dCA3dt = (CA2-CA3)/tau - k*CA3; dydt = [dCA1dt; dCA2dt; dCA3dt];
非线性微分方程边值问题打靶算法Matlab程序
![非线性微分方程边值问题打靶算法Matlab程序](https://img.taocdn.com/s3/m/b0efa93228ea81c759f57857.png)
非线性微分方程边值问题打靶算法Matlab程序【非线性微分方程边值问题打靶算法】参见:// matlabsky /thread-828-1-1.html【线性微分方程边值问题打靶算法】参见 :// matlabsky /thread-827-1-1.html 【线性微分方程边值问题有限差分算法】参见:// matlabsky /thread-829-1-1.html下面我们讲解下【非线性微分方程边值问题打靶算法】对于边值问题线性边值问题时,p、q和r只是x的函数,但是非线性时它们是x,y和y'的函数我们将问题转换为如下初值问题现在我们需要做的就是找到那个m,使得y(b)=beta,换句话说y现在由两个变量最终控制,即m和x。
对于这个求m的问题,我们可以使用牛顿迭代法得到(1)先给出一个初值m=1(2)计算出对应的y(b)为yb,这个直接使用ode45计算,得到y(end,1)就是yb(3)根据beta和yb的差异更新m(4)将新的m带入(2)重新计算yb并更新m(5)如此迭代直到m稳定或者在允许误差范围内下面详细描述了Matlab代码的运行过程以下内容需要回复才能看到复制内容到剪贴板代码:function [t,x]=nlineshoot(funcn,funcv,tspan,bound,tol,varargin)%对微分方程y''=p*y'+q*y+r,a<t<b,u(a)=alpha,u(b)=beta,其中p,q和r为x,y和y'的函数%只要对应修改p,q,r并将y'用x(2)替换,y用x(1)替换,就可以了%funcn=@(t,x)[x(2);p*x(2)+q*x(1)+r];%funcv=@(t,x)[x(2);p*x(2)+q*x(1)+r;x(4);x(3)*dF/dy+x(4)*dF/y'];%%% Example%F=y''=2y*y',y(0)=-1,y(pi/2)=1%%由方程有p=2y=2*x(1),q=r=0%dF/dy=2*y'=2*x(2) %y'用x(2)替换,y用x(1)替换%dF/dy'=2*y=2*x(1)%%funcn=@(t,x)[x(2);2*x(1)*x(2)];%funcv=@(t,x)[x(2);2*x(1)*x(2);x(4);2*x(2)*x(3)+2*x(1)*x(4)]; %tspan=[0 pi/2];%bound=[-1 1];%tol=1e-8;%[t,y]=nlineshoot(funcn,funcv,tspan,bound,tol);%plot(t,y)%legend('x1','x2')%%by dynamic%see also :// matlabsky%2009.3.12%%lb=bound(1);ub=bound(2);m0=0;m=1;while norm(m-m0)>tolm0=m;[t,x]=ode45(funcv,tspan,[lb;m;0;1],varargin);m=m0-(x(end,1)-ub)/x(end,3);end[t,x]=ode45(funcn,tspan,[lb;m]);下面的图形是代码中附带实例的运行结果线性微分方程边值问题打靶算法Matlab程序【非线性微分方程边值问题打靶算法】参见:// matlabsky /thread-828-1-1.html 【线性微分方程边值问题打靶算法】参见:// matlabsky /thread-827-1-1.html【线性微分方程边值问题有限差分算法】参见:// matlabsky /thread-829-1-1.html注意该算法只能完成二阶常微分方程双边值问题求解,至于其他形式的边值问题必须先转换到二阶形式对于下面的二阶常微分方程利用上面方程的线性结果和两个特殊的初值问题,我们可以构造两个等效的常微分初值方程初值问题,对于初值问题我们就可以直接使用ode**计算器或者龙哥库塔算法求解了。
Matlab求解边值问题方法+例题
![Matlab求解边值问题方法+例题](https://img.taocdn.com/s3/m/04008ee0b8f67c1cfad6b8f1.png)
y1 y , y2 y1
原方程组等价于以下标准形式的 方程组:
y y 1 2 y2 (1 x 2 ) y1 1
边界条件为:
y1 (0) 1 0 y1 (1) 3 0
边值问题的求解
求解:
c=1; solinit=bvpinit(linspace(0,4,10),[1 1]); sol=bvp4c(@ODEfun,@BCfun,solinit,[],c); x=[0:0.1:4]; y=deval(sol,x); plot(x,y(1,:),'b-',sol.x,sol.y(1,:),'ro') legend('解曲线','初始网格点解') % 定义ODEfun函数 function dydx=ODEfun(x,y,c) dydx=[y(2);-c*abs(y(1))]; % 定义BCfun函数 function bc=BCfun(ya,yb,c) bc=[ya(1);yb(1)+2];
z(0) 1, z (0) 0, z ( ) 0
本例中,微分方程与参数λ的数值有关。一般而言,对于任意的λ值,该问题无解, 但对于特殊的λ值(特征值),它存在一个解,这也称为微分方程的特征值问题。 对于此问题,可在bvpinit中提供参数的猜测值,然后重复求解BVP得到所需的 参数,返回参数为sol.parameters
边值问题的求解
如果取1,计算结果如何?
infinity=6; solinit=bvpinit(linspace(0,infinity,5),[0 0 1]); options=bvpset('stats','on'); sol=bvp4c(@ODEfun,@BCfun,solinit,options); eta=sol.x; f=sol.y; fprintf('\n'); 2 求解: f ff (1 ( f ) ) 0 fprintf('Cebeci & Keller report f''''(0)=0.92768.\n') fprintf('Value computed here is f''''(0)=%7.5f.\n',f(3,1)) 0.5 clf reset plot(eta,f(2,:)); axis([0 infinity 0 1.4]); title... f (0) 0, f (0) 0, ('Falkner-Skan equation, positive wall shear,\beta=0.5.') 边界条件: f ( ) 1, xlabel('\eta'),ylabel('df/d\eta'),shg % 定义ODEfun函数 function dfdeta=ODEfun(eta,f) beta=0.5; dfdeta=[f(2);f(3);-f(1)*f(3)-beta*(1-f(2)^2)]; % 定义BCfun函数 function res=BCfun(f0,finf) res=[f0(1);f0(2);finf(2)-1];
matlab求边值问题例题
![matlab求边值问题例题](https://img.taocdn.com/s3/m/857eab20a88271fe910ef12d2af90242a895abec.png)
matlab求边值问题例题摘要:一、引言二、MATLAB求边值问题的基本步骤1.建立微分方程模型2.离散化3.求解离散化后的线性方程组4.反演法求解边值问题三、实例分析1.二维热传导方程2.一维波动方程四、总结与展望正文:一、引言边值问题广泛应用于物理、工程、数学等领域,解决边值问题的一种有效方法是利用MATLAB进行数值求解。
本文将简要介绍MATLAB求边值问题的基本步骤,并通过实例分析二维热传导方程和一维波动方程的求解过程,以展示MATLAB在边值问题求解中的应用。
二、MATLAB求边值问题的基本步骤1.建立微分方程模型根据实际问题,建立相应的微分方程模型。
例如,对于二维热传导方程,可表示为:u/x = α * u/t2.离散化将空间和时间进行离散化,采用有限差分法将微分方程转化为线性方程组。
3.求解离散化后的线性方程组利用MATLAB内置的求解线性方程组的函数,如“ solve() ”,求解离散化后的线性方程组。
4.反演法求解边值问题对于具有周期性的边值问题,可以采用反演法进行求解。
具体步骤如下:1) 求解对应的周期边界条件下的一维线性方程组;2) 利用MATLAB的插值函数,如“interp1()”或“spline()”,对解进行插值;3) 通过反演公式,求解原始边值问题。
三、实例分析1.二维热传导方程考虑如下二维热传导方程的边值问题:u/x = α * u/t,u(x, 0) = f(x),u(0, t) = g(t),u(a, t) = u(b, t)根据基本步骤,首先建立微分方程模型,然后进行离散化,利用MATLAB 求解离散化后的线性方程组。
最后,通过反演法求解边值问题。
2.一维波动方程考虑如下边值问题:i * u/t = u/x,u(x, 0) = f(x),u(a, t) = u(b, t)同样地,根据基本步骤,建立微分方程模型,离散化,求解离散化后的线性方程组,最后采用反演法求解边值问题。
(偏)微分方程一类边值问题的数值求解(附matlab程序)
![(偏)微分方程一类边值问题的数值求解(附matlab程序)](https://img.taocdn.com/s3/m/39326bd849649b6648d747c5.png)
(偏)微分方程一类边值问题的数值求解本文介绍了椭圆型(偏)微分方程一类边值问题的数值求解程序(笔者自编)及其使用方法。
程序基于差分原理,将连续的(偏)微分方程在节点上离散化,最终化为线性代数方程组。
对于原理方面的问题很多微分方程数值解的参考书中都有详尽的描述,但一般缺少具体的实现程序。
笔者认为,对这类数值方法是否理解并能运用的关键之一还是在于是否能写出用于解决问题的程序。
理论基础固然必不可少,程序确往往是我们解决问题的敲门砖。
一维椭圆型方程可表示如下:],[,0)(,0)()(b a x x q x p f qu dxdur dx du p dx d Lu ∈>>=++-=其中L 表示微分算子,很明显这是一个线性算子。
这里要求p ,q 大于零是为了保证最终得到的线性方程组有唯一的非零解,但事实上不满足这个条件可能也是有解的,这涉及到微分方程解的存在性也确定性问题,读者若有兴趣可参考相关书籍。
更具体的,我们可以举出一个一维椭圆型方程的例子来: 例(1):)cos(10)exp(222x u x dx dux dxu d =++- 即:)cos(10),ex p()(,)(,1)(2x f x x q x x r x p ====⎩⎨⎧==3)3(1)1(u u 利用文中提及的程序,可以将以上问题表述为:syms x %定义符号变量xp=@(x) 1; r=@(x) x.^2; q=@(x) exp(x); f=@(x) 10*cos(x);Odbvp(p,r,q,f,1,3,1,3);11.2 1.4 1.6 1.82 2.2 2.4 2.6 2.83-0.50.511.522.532~Odbvp u=f(x)xu哈哈,一条非常漂亮的曲线。
若不满足p,q>0,我们也举一例: 例(2)syms xp=@(x) 10*cos(x); r=@(x) x.^2;q=@(x) 10*sin(x); f=@(x) 10*cos(x);Odbvp(p,r,q,f,1,3,1,3)1 1.2 1.4 1.6 1.82 2.2 2.4 2.6 2.83-500501001502002503003502~Odbvp u=f(x)xu可以发现,程序仍能求解,但结果的光滑性不好。
偏微分方程边值问题的数值解法
![偏微分方程边值问题的数值解法](https://img.taocdn.com/s3/m/749669fcc8d376eeaeaa31a1.png)
求解偏微分方程的边值问题本实验学习使用MATLAB 的图形用户命令pdetool 来求解偏微分方程的边值问题。
这个工具是用有限元方法来求解的,而且采用三角元。
我们用内个例题来说明它的用法。
一、MATLAB 支持的偏微分方程类型考虑平面有界区域D 上的二阶椭圆型PDE 边值问题:()c u u f α-∇∇+= (1.1)其中 (1) , (2) a,f D c x y ⎛⎫∂∂∇=⨯ ⎪∂∂⎝⎭是上的已知函数(3)是标量或22的函数方阵未知函数为(,) (,)u x y x y D ∈。
它的边界条件分为三类:(1)Direchlet 条件:hu f = (1.2)(2)Neumann 条件: ()n c u qu g ∇+= (1.3)(3)混合边界条件:在边界D ∂上部分为Direchlet 条件,另外部分为Neumann 条件。
其中,,,,h r q g c 是定义在边界D ∂的已知函数,另外c 也可以是一个2*2的函数矩阵,n 是沿边界的外法线的单位向量。
在使用pdetool 时要向它提供这些已知参数。
二、例题例题1 用pdetool 求解 22D 1 D: 10u x y u ∂⎧-∆=+≤⎪⎨=⎪⎩ (1.4)解 :首先在MATLAB 的工作命令行中键入pdetool ,按回牟键确定,于是出现PDE Toolbox 窗口,选Genenic Scalar 模式.( l )画区域圆单击椭圆工具按钮,大致在(0,0)位置单击鼠标右键,拖拉鼠标到适当位置松开。
为了保证所绘制的圆是标准的单位园,在所绘园上双击,打开 Object Dialog 对话框,精确地输入圆心坐标X-center 为0 、Y-center 为0 及半径Radius 为l ,然后单击OK 按钮,这样单位画已画好.( 2 )设置边界条件单击工具边界模式按钮 ,图形边界变红,逐段双击边界,打开Boundary condition 对话框.输入边界条件.对于同一类型的边界,可以按Shift 键,将多个边界同时选择,统一设边界条件.本题选择Dirichlet 条件,输入h 为1 , r 为0。
matlab求解初边值问题的偏微分方程
![matlab求解初边值问题的偏微分方程](https://img.taocdn.com/s3/m/da972d9dac51f01dc281e53a580216fc710a5364.png)
偏微分方程是描述自然界中动态过程的重要数学工具,在工程领域中,求解偏微分方程是很多实际问题的重要一步。
MATLAB作为一种强大的数学计算软件,提供了丰富的工具和函数来求解各种类型的偏微分方程。
本文将介绍使用MATLAB求解初边值问题的偏微分方程的方法和步骤。
一、MATLAB中的偏微分方程求解工具MATLAB提供了几种可以用来求解偏微分方程的工具和函数,主要包括:1. pdepe函数:用于求解偏微分方程初边值问题的函数,可以处理各种类型的偏微分方程,并且可以自定义边界条件和初始条件。
2. pdepeplot函数:用于绘制pdepe函数求解得到的偏微分方程的解的可视化图形,有助于直观地理解方程的解的特性。
3. pdetool工具箱:提供了一个交互式的图形用户界面,可以用来建立偏微分方程模型并进行求解,适用于一些复杂的偏微分方程求解问题。
二、使用pdepe函数求解偏微分方程初边值问题的步骤对于给定的偏微分方程初边值问题,可以按照以下步骤使用pdepe函数进行求解:1. 定义偏微分方程需要将给定的偏微分方程转化为标准形式,即将偏微分方程化为形式为c(x,t,u)∂u/∂t = x(r ∂u/∂x) + ∂(p∂u/∂x) + f(x,t,u)的形式。
2. 编写边界条件和初始条件函数根据实际问题的边界条件和初始条件,编写相应的函数来描述这些条件。
3. 设置空间网格选择合适的空间网格来离散空间变量,可以使用linspace函数来产生均匀分布的网格。
4. 调用pdepe函数求解偏微分方程将定义好的偏微分方程、边界条件和初始条件函数以及空间网格作为参数传递给pdepe函数,调用该函数求解偏微分方程。
5. 可视化结果使用pdepeplot函数绘制偏微分方程的解的可视化图形,以便对解的性质进行分析和理解。
三、实例分析考虑一维热传导方程初边值问题:∂u/∂t = ∂^2u/∂x^2, 0<x<1, 0<t<1u(0,t) = 0, u(1,t) = 0, u(x,0) = sin(πx)使用MATLAB求解该初边值问题的步骤如下:1. 定义偏微分方程将热传导方程化为标准形式,得到c(x,t,u) = 1, r = 1, p = 1, f(x,t,u) = 0。
matlab求边值问题例题
![matlab求边值问题例题](https://img.taocdn.com/s3/m/d228f9acdbef5ef7ba0d4a7302768e9950e76e4b.png)
matlab求边值问题例题
【最新版】
目录
1.MATLAB 求边值问题的基本概念
2.边值问题的例题解析
3.MATLAB 求解边值问题的步骤
4.总结
正文
一、MATLAB 求边值问题的基本概念
边值问题是数学物理中常见的一种问题,它是指在给定的边界条件下,求解偏微分方程的数值解。
在 MATLAB 中,我们可以通过边界元方法或者有限元方法来求解这类问题。
二、边值问题的例题解析
假设我们有一个二维的 Laplace 方程,给定边界条件为:x=0 时,y"=0;x=1 时,y"=0;y=0 时,x"=0;y=1 时,x"=0。
我们可以通过以下步骤在 MATLAB 中求解该问题。
1.定义边界和区域:在 MATLAB 中,我们首先需要定义问题的边界和区域。
我们可以使用边界元方法,将边界和区域定义为一个二维矩阵。
2.定义方程:然后,我们需要定义 Laplace 方程。
在 MATLAB 中,我们可以使用 PDE 工具箱中的函数来定义方程。
3.定义边界条件:接下来,我们需要定义边界条件。
在 MATLAB 中,我们可以使用边界元方法中的函数来定义边界条件。
4.求解:最后,我们可以使用 MATLAB 中的函数来求解该问题。
我们可以使用 PDE 工具箱中的函数,或者使用我们自己编写的函数来求解。
三、MATLAB 求解边值问题的步骤
1.定义边界和区域
2.定义方程
3.定义边界条件
4.求解
四、总结
MATLAB 是一个强大的数学软件,它可以帮助我们求解各种各样的数学问题,包括边值问题。
数学实验“微分方程组边值问题数值算法(打靶法,有限差分法)”实验报告(内含matlab程序)
![数学实验“微分方程组边值问题数值算法(打靶法,有限差分法)”实验报告(内含matlab程序)](https://img.taocdn.com/s3/m/8d9eb78e50e79b89680203d8ce2f0066f53364ae.png)
西京学数学软件实验任务书课程名称数学软件实验班级数0901学号0912020107姓名李亚强微分方程组边值问题数值算法(打靶法,有限差分法)实验课题熟悉微分方程组边值问题数值算法(打靶法,有限差实验目的分法)运用Matlab/C/C++/Java/Maple/Mathematica等其中实验要求一种语言完成微分方程组边值问题数值算法(打靶法,有限差分法)实验内容成绩教师实验二十七实验报告1、实验名称:微分方程组边值问题数值算法(打靶法,有限差分法)。
2、实验目的:进一步熟悉微分方程组边值问题数值算法(打靶法,有限差分法)。
3、实验要求:运用Matlab/C/C++/Java/Maple/Mathematica 等其中一种语言完成程序设计。
4、实验原理:1.打靶法:对于线性边值问题(1)⎩⎨⎧==∈=+'+''βα)(,)(],[)()()(b y a y b a x x f y x q y x p y 假设是一个微分算子使:L ()()Ly y p x y q x y '''=++则可得到两个微分方程:,,)(1x f Ly =α=)(1a y 0)(1='a y ,, (2)⇔)()()(111x f y x q y x p y =+'+''α=)(1a y 0)(1='a y ,,02=Ly 0)(2=a y 1)(2='a y ,, (3)⇔0)()(222=+'+''y x q y x p y 0)(2=a y 1)(2='a y 方程(2),(3)是两个二阶初值问题.假设是问题1y(2)的解,是问题(3)的解,且,则线性边值问2y 2()0y b ≠题(1)的解为: 。
1122()()()()()y b y x y x y x y b β-=+2.有限差分法:基本思想是把连续的定解区域用有限个离散点构成的网格来代替,这些离散点称作网格的节点;把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似;把原方程和定解条件中的微商用差商来近似, 积分用积分和来近似,于是原微分方程和定解条件就近似地代之以代数方程组,即有限差分方程组 , 解此方程组就可以得到原问题在离散点上的近似解。
边值问题(BVP)的Matlab解法
![边值问题(BVP)的Matlab解法](https://img.taocdn.com/s3/m/393065394a73f242336c1eb91a37f111f1850d0f.png)
边值问题(BVP)的Matlab解法求微分⽅程(1+x)D2y=2y-4 初始条件 y(0)=0 y(1)=2Dy(1)如果想⽤inline和ode45解决,不⽤function窗⼝,该如何做?解1. 对于此类边值微分⽅程,ode**函数是⽆⼒直接求解的,Matlab提供了bvp解算器。
2. 对⼲你⽤的Matlab版本,⽤@(x)函数(匿名函数)⽐inline更⽅便。
3. 请参考本⽂中的其它例题及相关资料理解下⾯的代码。
4. 在Matlab7.1版上,可⽤: dsolve('D2y = (2*y-4)/(1+t)', 'y(0) = 0', 'Dy(1) = y(1)/2') 求出解析解(符号解)。
如不需要解析解与数值解的对⽐,可不要第⼀段代码。
% 解析解 syms t y = -(1+t)^(1/2)*besseli(1,2*2^(1/2)*(1+t)^(1/2))*(4*i*bessely(0,4*i)-2*bessely(1,4*i)+2^(1/2)*bessely(1,2*i*2^(1/2)))/(2*i*bessely(0,4*i)*besseli(1,2*2^(1/2))-bessely(1,4*i)*besseli(1,2*2^(1/2))+besseli(1,4)*bessely(1,2*i*2^(1/2))-2*besseli(0,4)*bessely(1,2*i*2^(1/2)))+(1+t)^(1/2)*bessely(1,2*i*2^(1/2)*(1+t)^(1/2))*(besseli(1,2*2^(1/2))*2^(1/2)-2*besseli(1,4)+4*besseli(0,4))/(2*i*bessely(0,4*i)*besseli(1,2*2^(1/2))-bessely(1,4*i)*besseli(1,2*2^(1/2))+besseli(1,4)*bessely(1,2*i*2^(1/2))-2*besseli(0,4)*bessely(1,2*i*2^(1/2)))+2; ezplot(y,[0 1])% 数值解 dydx = @(x,y) [y(2);(2*y(1)-4)/(1+x) ]; %边值微分⽅程 res = @(ya,yb) [yb(2) - yb(1)/2;ya(1) - 0 ]; %边界条件 solinit =bvpinit(linspace(0,1,10),[1 0]); sol = bvp4c(dydx,res,solinit); xint = linspace(0,1,50); Sxint = deval_r(sol,xint);% 画图 hold on plot(xint,Sxint(1,:),'*') title('stevenchang041''s equation.') xlabel('x') ylabel('solution y')。
两点边值问题的不同迭代法比较及matlab实现
![两点边值问题的不同迭代法比较及matlab实现](https://img.taocdn.com/s3/m/03e6244765ce05087732138a.png)
两点边值问题的不同迭代法比较及matlab实现问题:考虑两点边值问题:2,dydyaa,,,0,,1,,2 dxdx,,yy(0),0,(1),1.,x,1,a,容易知道它的精确解为: y,(1,e),ax1/,,1,ex,ih,i,1,2,...n,1为了将微分方程离散,把[0,1]区间n等分,令h=1/n,,得到差分方程 i2(,,h)y,(2,,h)y,,y,ah,从而得到迭代方程组的系数矩阵A。
i,ii,11 对=1,a=1/2,n=100,分别用jacobi,G-S,超松弛迭代法分别求线性方程组的解,要求,4位有效数字,然后比较与精确解的误差。
对=0.1,=0.01,=0.001,考虑同样问题。
,,,思想:利用书上的迭代公式即可。
注意问题:迭代矩阵是n-1阶的,不是n阶;等号右端向量b的最后一项,不是ah^2,而是ah^2-eps-hx,1,a,精确解: y,(1,e),ax1/,,1,e带入a=1/2,=1 ,代码:>> clear>> x=linspace(0,1); truy=(1-0.5)/(1-exp(-1/1))*(1-exp(-x./1))+x.*0.5;figure;plot(x,truy,'g','LineWidth',1.5); hold on;Grid图:三种方法的实现Jacobi法:代码见附录Eps=1结果:迭代次数k:22273结果与精确解的比较图(绿色粗线是精确解,黑色细线是迭代结果)Eps=0.1结果:迭代次数k:8753结果与精确解的比较图(绿色粗线是精确解,黑色细线是迭代结果) Eps=0.01结果:迭代次数k:661结果与精确解的比较图(绿色粗线是精确解,黑色细线是迭代结果)G-S迭代法:代码见附录Eps=1结果:迭代次数k:11125结果与精确解的比较图(绿色粗线是精确解,黑色细线是迭代结果)Eps=0.1结果:迭代次数k:4394结果与精确解的比较图(绿色粗线是精确解,黑色细线是迭代结果)Eps=0.01结果:迭代次数k:379结果与精确解的比较图(绿色粗线是精确解,黑色细线是迭代结果)超松弛法:代码见附录Eps=1 w=1.56结果:迭代次数k:3503结果与精确解的比较图(绿色粗线是精确解,黑色细线是迭代结果)Eps=0.1 w=1.56结果:迭代次数k:1369结果与精确解的比较图(绿色粗线是精确解,黑色细线是迭代结果)Eps=0.01 w=1.56结果:迭代次数k:131结果与精确解的比较图(绿色粗线是精确解,黑色细线是迭代结果) 结果:Jacobi、G-S、超松弛法,三者都能够取得对精确解的良好逼近,但是,在相同的精度条件下,三者的收敛速度是不一样的,jacobi<G-S<超松弛,也就是说,在迭代次数相同的条件下,精度:jacobi<G-S<超松弛。
李荣华二维抛物方程的初边值问题matlab
![李荣华二维抛物方程的初边值问题matlab](https://img.taocdn.com/s3/m/9e16452b24c52cc58bd63186bceb19e8b8f6ecff.png)
李荣华二维抛物方程的初边值问题在数学和工程领域中扮演着重要角色。
通过使用Matlab,我们能够更深入地理解和解决这一问题。
在本文中,我们将从基础概念开始,逐步深入讨论李荣华二维抛物方程的初边值问题,并结合Matlab进行实际分析和解决。
一、初边值问题的概念和涵义1. 初边值问题的定义初边值问题是指在偏微分方程中,除了部分边界条件外,还需给出一些初始条件。
李荣华二维抛物方程的初边值问题即是在方程中给定初值条件和边界条件的问题。
2. 李荣华二维抛物方程简介李荣华二维抛物方程是描述热传导、扩散等现象的数学模型,在物理学和工程领域有着广泛的应用。
它的形式通常为一个关于未知函数和其在空间和时间上的导数的方程。
3. Matlab在初边值问题中的应用Matlab作为一种强大的数学建模和仿真工具,能够帮助我们求解各种类型的偏微分方程,包括初边值问题。
通过Matlab,我们可以对李荣华二维抛物方程进行数值求解和模拟,从而更清晰地理解问题本质。
二、李荣华二维抛物方程的初边值问题解析1. 方程的建立我们需要建立李荣华二维抛物方程的数学模型,并给定相应的初值条件和边界条件。
在Matlab中,我们可以通过定义矩阵和方程的形式来实现这一过程。
2. 数值求解我们利用Matlab提供的数值求解方法,如有限差分法或有限元法,对初边值问题进行求解。
通过程序的编写和运行,我们可以得到方程的数值解,并对物理过程有更直观的认识。
3. 结果分析在求得数值解之后,我们可以对结果进行分析和可视化。
通过Matlab 的绘图功能,我们能够直观地观察李荣华二维抛物方程的演化过程,以及初值条件和边界条件对解的影响。
三、个人观点和总结在本文中,我们深入探讨了李荣华二维抛物方程的初边值问题,并结合Matlab进行了实际分析和解决。
通过对问题的全面评估和数值求解,我们对方程的特性和行为有了更深入的理解。
个人观点上,我认为Matlab是一个非常强大的工具,能够帮助我们更直观、高效地处理数学问题。
二阶椭圆偏微分方程实例求解(附matlab代码)
![二阶椭圆偏微分方程实例求解(附matlab代码)](https://img.taocdn.com/s3/m/5aa10dd3b9d528ea80c77947.png)
《微分方程数值解法》期中作业实验报告二阶椭圆偏微分方程第一边值问题姓名:学号:班级:2013年11月19日二阶椭圆偏微分方程第一边值问题摘要对于解二阶椭圆偏微分方程第一边值问题.课本上已经给出了相应的差分方程。
而留给我的难题就是把差分方程组表示成系数矩阵的形式.以及对系数进行赋值。
解决完这个问题之后.我在利用matlab 解线性方程组时.又出现“out of memory ”的问题。
因为99*99阶的矩阵太大.超出了分配给matlab 的使用内存。
退而求其次.当n=10.h=1/10或n=70.h=1/70时.我都得出了很好的计算结果。
然而在解线性方程组时.无论是LU 分解法或高斯消去法.还是gauseidel 迭代法.都能达到很高的精度。
关键字:二阶椭圆偏微分方程 差分方程 out of memory LU 分解 高斯消去法 gauseidel 迭代法一、题目重述解微分方程:()()2222((,))((,))()(,)()(,)(,)1y x x x y y x y yxxyxye u x y e u x y x y u x y x y u x y u x y y e x e e y x e--+++-+=-++++已知边界:(0,)1,(1,),(,0)1,(,1)y x u y u y e u x u x e ====求数值解, 把区域[0,1][0,1]G =?分成121/100,1/100h h ==.n =100 注:老师你给的题F 好像写错了.应该把22x y y e x e +改成22y x y e x e +。
二、问题分析与模型建立2.1微分方程上的符号说明()()22221y x xy xy y e x e e y x e -++++2.2课本上差分方程的缺陷课本上的差分方程为:举一个例子:当i=2,j=3时.;当i=3,j=3时.。
但是.显然这两个不是同一个数.其大小也不相等。
一类基于MATLAB的微分方程边值问题数值解的算法研究
![一类基于MATLAB的微分方程边值问题数值解的算法研究](https://img.taocdn.com/s3/m/662f37f3541810a6f524ccbff121dd36a32dc46b.png)
一类基于MATLAB的微分方程边值问题数值解的算法研究王来英;薛亚宏
【期刊名称】《中国西部科技》
【年(卷),期】2012(011)008
【摘要】本文对线性与非线性微分方程边值问题的打靶算法进行了理论分析,进而建立了解决一般微分方程边值问题的数学模型,并通过MATLAB来实现微分方程中的一类边值问题的计算机求解,使其具备更好的可操作性与变量的可扩展性.同时利用MATLAB强大的图形绘制功能提供良好的可视化环境;最后以实例形式对MATLAB算法进行实变量运行与观测.
【总页数】3页(P48-50)
【作者】王来英;薛亚宏
【作者单位】甘肃省天水市职业技术学校数学教研室,甘肃天水741000;甘肃工业职业技术学院电信学院,甘肃天水741025
【正文语种】中文
【相关文献】
1.两类基于MATLAB的非线性微分方程数值解的算法研究 [J], 薛亚宏
2.基于Matlab常微分方程数值解的分析与比较 [J], 马泽涛;文鹏;施琳达;施玲玲;刘耀桓;吝维军
3.一类脉冲时滞偏微分方程初边值问题的matlab图像 [J], 邹序焱;饶若峰
4.两类基于MATLAB的非线性微分方程数值解的算法研究 [J], 薛亚宏
5.一类二阶线性常微分方程两点边值问题的数值解 [J], 雍龙泉;刘三阳;熊文涛;迟晓妮
因版权原因,仅展示原文概要,查看原文内容请购买。
一类基于MATLAB的微分方程边值问题数值解的算法研究_王来英
![一类基于MATLAB的微分方程边值问题数值解的算法研究_王来英](https://img.taocdn.com/s3/m/3070b5aeb0717fd5360cdc9b.png)
3
微分方程边值问题的打靶算法
3.1 线性微分方程边值问题的打靶算法 对于线性微分方程: ( x) q( x) y( x) f ( x) , y ( x) p( x) y 其中 p(x)、q(x)和 f(x)均为给定函数,则上面的边值条件可简 化为:y(a)=ra,y(b)=rb。 而打靶算法(shooting method)的思想是找到能够满足 然后再利用常规初值算 边界点条件的相应初值 y(0)和 y (0) , 法去求解这个初值问题[2]。其由以下五个步骤来完成: 首先,求出方程初值问题的数值解 y1(b)、y2(b)、yp(b), 其中若 y2(b)≠0,则计算: m b a y1 (b) y p (b) 。
(e2 -3) e x +(3-e) e2 x 3 1 + + x 4 e(e-1) 4 2
其计算的精度可由下面语句得出:>>y0=((exp(2)-3)* exp(t)+ (3-exp(1)*exp(2*t))/(4*exp(1)-1))+3/4+t/2;e1=norm
中国西部科技
2012 年 08 月第 11 卷第 08 期总第 277 期
y 2 (b )
图1 线性微分方程边值问题图形解
事实上由微分方程理论可得原方程的解析解为:
y(x) =
然后计算下面初值问题的数值解,则 y ( x) 即为原边值
( x) p( x) y ( x) q ( x) y ( x) f ( x) , y ( a ) ra , 问题的数值解: y (a) m 。 y
a a , b y (b) b y (b) b ,则需要给出更一般的 a y
matlab三边测量算法横纵坐标图解=
![matlab三边测量算法横纵坐标图解=](https://img.taocdn.com/s3/m/06ab3837a5e9856a561260a1.png)
球面短程线计算公式
L R
①输入经纬度数据和地球半径; ②转换两城市的经纬度为地心直角坐标数据; ③提取两个点的向径坐标; ④计算向径间的夹角和短程线长度并输出计算结果。
10/22
按顺序录入程序文件(文件名:distance.m)
16/22
————程序设计中的流程控制————
例2.14 3n + 1 问题. 对任一自然数n,按如下法则进行运算:若n为偶数, 则将n除2;若n为奇数,则将n乘3加1。将运算结果按 上面法则继续运算, 重复若干次后结果最终是1.
n=input(‘input n=’); %输入数据 while n~=1 r=rem(n,2); %求n/2的余数 if r = =0 n=n/2 %第一种操作 else n=3*n+1 %第二种操作 end end n=5 16, 8, 4, 2, 1
7/22
例2.3直线平行于Z轴沿x-y平面上的四边形移动,形 成四边形柱面。利用矩阵方法绘制四边形柱面.
设四边形顶点为(-1, -1), (1, -1), (1, 1), (-1, 1).由于四边 形是封闭的图形,将第五个点设为第一个点。设柱面 高为1,创建Z坐标矩阵
0 0 0 0 0 Z 1 1 1 1 1
R=6400; S0=4*pi*R*R; d=10000; S=2*pi*R*R*d/(R+d); S/S0*100
6/22
M文件分为命令文件和函数文件两种
命令文件——MATLAB命令的有序集合。 文件执行——对文件中命令进行批处理,即从第一 条命令开始按顺序执行,直到最后一条命令。如果 中间某条命令出错,则中断并输出错误信息 ①在编辑窗口中编写; ②保存并对文件命名; ③命令窗口键入文件 名运行; ④观察运行结果; · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
matlab有限元分析实例
![matlab有限元分析实例](https://img.taocdn.com/s3/m/0471846a6137ee06eef9182e.png)
matlab有限元分析实例问题描述:考虑一平面有界区域,设其边界为[。
我们求解泊松方程之狄利克雷边值问题。
问题的强形式为一椭圆型偏微分方程当之几何形状稍复杂时,一般无法求得其解析解。
我们可应用有限元法来求其数值解。
通常我们先写出该问题的抽象弱形式:求使得其中为检验函数为一适当索伯列夫空间(在本例中也是一希尔伯特空间)为一双线性型为一线性型。
其具体表达式为有限元空间离散我们采用最简单的二维单元离散单元即三节点线性三角形单元,其插值基函数(即形函数为一次多项式。
有限元之核心思想为:使用离散的函数空间来分片逼近连续的函数空间。
于是所求之近似解可写成基函数之线性组合其中为待求系数,常称为自由度。
将此近似解之表达式代入前述问题之弱形式,并取检验函数,可得写成矩阵形式,便成为常见之有限元方程由于历史原因,通常采用固体力学中的习惯命名:[公式] 为刚度矩阵,[公式] 为自由度向量,[公式] 为载荷向量。
由于空间已分片离散,上面的有限元方程只在各单元内部成立。
为了求解方便,通常我们将所有单元的有限元方程连立起来求解,于是需要将各单元之刚度矩阵组装成总体刚度矩阵。
求解器MATLAB 代码及简释本求解器之十行代码如下:function u=fem(nds,els,bcs)nnd=size(nds,1); u=zeros(nnd,1); K=zeros(nnd,nnd); f=zeros(nnd,1);for j=1:size(els,1)K(els(j,:),els(j,:))=K(els(j,:),els(j,:))+stima(nds(els(j,:),:));f(els(j,:))=f(els(j,:))+ones(3,1)*det([1,1,1;nds(els(j,:),:)'])/ 6;endfreends=setdiff(1:nnd,bcs);u(freends)=K(freends,freends)\f(freends);function stima=stima(vertices)ndim=size(vertices,2);J=[ones(1,ndim+1);vertices'];B=J\[zeros(1,ndim);eye(ndim)];stima=det(J).*B*B'/prod(1:ndim);输入数据格式:唯一的一个函数stima 用来计算各单元刚度矩阵。
matlab用差分法求解边值问题
![matlab用差分法求解边值问题](https://img.taocdn.com/s3/m/0b875191b8f3f90f76c66137ee06eff9aef8499d.png)
MATLAB是一种强大的数学软件工具,用于解决各种数学问题,包括边值问题。
差分法是一种常用的数值分析方法,可以用于求解微分方程的边值问题。
本文将介绍如何使用MATLAB的差分法求解边值问题,包括算法原理、具体步骤和示例代码。
1. 算法原理差分法是一种数值解微分方程的常用方法,它的基本思想是用离散点代替连续函数,然后利用差分代替导数,将微分方程转化为代数方程,最终得到离散点上的数值解。
对于边值问题,差分法可以通过离散化区域,并在区域边界处使用相应的差分格式来求解。
2. 具体步骤(1)离散化区域需要将求解的区域进行离散化,将其划分为若干个离散点,形成一个网格。
这可以通过在区域内部建立均匀或非均匀的网格来实现。
(2)建立代数方程利用差分公式将原微分方程转化为代数方程。
根据边值问题的具体条件,可以构建相应的代数方程组。
(3)求解代数方程组利用MATLAB的线性方程求解器,如“backslash”运算符(\),求解得到代数方程组的解,即为边值问题的数值解。
3. 示例代码下面是一个简单的边值问题求解的MATLAB示例代码,以二阶常微分方程为例:``` MATLAB定义区域范围和离散步长a = 0;b = 1;N = 10; 离散点个数h = (b - a) / N; 离散步长构建代数方程组A = zeros(N-1, N-1);b = zeros(N-1, 1);for i = 1:N-1A(i, i) = 2/h^2;if i > 1A(i, i-1) = -1/h^2;endif i < N-1A(i, i+1) = -1/h^2;endb(i) = f(a + i*h); 右端项end求解代数方程组u = A \ b;绘制数值解x = a:h:b;plot(x, [0; u; 0]);```在这个示例代码中,我们首先定义了求解区域范围和离散步长,然后根据差分格式建立了代数方程组,并使用MATLAB的线性方程求解器求解得到数值解,最后绘制了数值解的图像。
matlab求边值问题例题
![matlab求边值问题例题](https://img.taocdn.com/s3/m/cc411679590216fc700abb68a98271fe900eaf55.png)
matlab求边值问题例题
一个常见的边值问题是求解方程 -d^2u/dx^2 = f(x),其中
d^2u/dx^2是二阶导数,f(x)是已知函数,边值条件是u(0) = a 和u(L) = b,L是给定常数。
这个方程可以用来描述一个细杆在两端受到固定边界条件的情况下的温度分布问题。
下面是一个解决这个边值问题的MATLAB例子:
步骤1:定义相关参数和函数
L = 1; % 简化计算,取L = 1
a = 0; % 边值条件 u(0) = 0
b = 2; % 边值条件 u(L) = 2
f = @(x) x; % 定义已知函数 f(x) = x
步骤2:离散化求解区域
n = 100; % 离散化步数
x = linspace(0, L, n); % 在[0, L]区间均匀离散化
步骤3:定义线性方程组
A = zeros(n, n); % 构建系数矩阵A
b = zeros(n, 1); % 构建右端项b
h = L/n; % 步长
% 内部节点的方程
for i = 2:n-1
A(i, i) = -2/h^2;
A(i, i-1) = 1/h^2;
A(i, i+1) = 1/h^2;
b(i) = f(x(i));
end
% 边界条件
A(1, 1) = 1;
A(n, n) = 1;
b(1) = a;
b(n) = b;
步骤4:解线性方程组u = A\b;
步骤5:绘制结果plot(x, u)
xlabel('x')
ylabel('u')
title('边值问题解')。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
边值问题求解指令:bvp4c()
输出参数: 输出变量sol是一个结构体 sol.x是指令bvp4c所采用的网格节点; sol.y是y(x)在sol.x网点上的近似解值; sol.yp是y'(x)在sol.x网点上的近似解值; sol.parameters是微分方程所包含的未知参数的近似解值。 当被解微分方程包含未知参数时,该域存在。
使残差不断减小,从而获得满足精度要求的解
Matlab求解边值问题的基本指令
solinit=bvpinit(x,v,parameters) 生成bvp4c调用指令所必须的“解猜测网” sol=bvp4c(odefun,bcfun,solinit,options,题的求解
y y 10 求解两点边值问题: y (0) y(1) 0
令:
y1 y , y2 y1
原方程组等价于以下标准形式的 方程组:
y1 y2 y2 y1 10
边界条件为:
y1 (0) 0, y1 (1) 0
边值问题求解指令:bvp4c()
sol=bvp4c(odefun,bcfun,solinit,options,p1,p2,…) 输入参数:
odefun是计算导数的m函数文件。该函数的基本形式为: dydx=odefun(x,y,parameters,p1,p2,…),在此,自变量x是标量,y, dydx是列向量。 bcfun是计算边界条件下残数的m函数文件。其基本形式为: res=bcfun(ya,yb,parameters,p1,p2,…),文件输入宗量ya,yb是边界 条件列向量,分别代表y在a和b处的值。res是边界条件满足时的残数 列向量。注意:例如odefun函数的输入宗量中包含若干“未知”和 “已知”参数,那么不管在边界条件计算中是否用到,它们都应作为 bcfun的输入宗量。 输入宗量options是用来改变bvp4c算法的控制参数的。在最基本用法 中,它可以缺省,此时一般可以获得比较满意的边值问题解。如需更 改可采用bvpset函数,使用方法同odeset函数。 输入宗量p1,p2等表示希望向被解微分方程传递的已知参数。如果无须 向微分方程传递参数,它们可以缺省。
边值问题的求解
y (1 x 2 ) y 1 求解: y (0) 1, y (1) 3
令:
solinit=bvpinit(linspace(0,1,10),[0 1]); sol=bvp4c(@ODEfun,@BCfun,solinit); x=[0:0.1:1]; y=deval(sol,x); xP=[0:0.1:1.0]; yP=[1 1.0743 1.1695 1.2869 1.4284... 1.5965 1.7947 2.0274 2.3004 2.6214 3]; plot(xP,yP,'o',x,y(1,:),'r-') legend('Analytical Solution','Numerical Solution',... 'location','Northwest') legend boxoff % 定义ODEfun函数 function dydx=ODEfun(x,y) dydx=[y(2);(1+x^2)*y(1)+1]; % 定义BCfun函数 function bc=BCfun(ya,yb) bc=[ya(1)-1;yb(1)-3];
solinit=bvpinit(linspace(0,1,10),[1 0]); sol=bvp4c(@ODEfun,@BCfun,solinit); x=[0:0.05:0.5]; y=deval(sol,x); xP=[0:0.1:0.5]; yP=[0 -0.41286057 -0.729740656... -0.95385538 -1.08743325 -1.13181116]; plot(xP,yP,'o',x,y(1,:),'r-') legend('Analytical Solution','Numerical Solution') % 定义ODEfun函数 function dydx=ODEfun(x,y) dydx=[y(2);y(1)+10]; % 定义BCfun函数 function bc=BCfun(ya,yb) bc=[ya(1);yb(1)];
sxint=deval(sol,xint) 计算微分方程积分区间内任何一点的解值
初始解生成函数:bvpinit()
solinit=bvpinit(x,v,parameters)
x指定边界区间[a,b]上的初始网络,通常是等距排列的(1×M)一维数组。 注意:使x(1)=a,x(end)=b;格点要单调排列。 v是对解的初始猜测 solinit(可以取别的任意名)是“解猜测网(Mesh)”。 它是一个结构体,带如下两个域: solinit.x是表示初始网格有序节点的(1×M)一维数组,并且 solinit.x(1)一定是a,solinit.x(end)一定是b。M不宜取得太大,10 数量级左右即可。 solinit.y是表示网点上微分方程解的猜测值的(N×M)二维数组。 solinit.y(:,i)表示节点solinit.x(i)处的解的猜测值。
Matlab求解边值问题方法:bvp4c函数
1. 把待解的问题转化为标准边值问题
y f ( x, y) g ( y(a), y(b)) 0
2. 因为边值问题可以多解,所以需要为期望解指定一个初始
猜测解。该猜测解网(Mesh)包括区间[a, b]内的一组网 点(Mesh points)和网点上的解S(x) 3. 根据原微分方程构造残差函数 r ( x) S ( x) f ( x, S ( x)) 4. 利用原微分方程和边界条件,借助迭代不断产生新的S(x),