Matlab求解边值问题方法+例题

合集下载

重要:MATLAB常微分方程(组)数值解法

重要:MATLAB常微分方程(组)数值解法

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程序

非线性微分方程边值问题打靶算法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求解边值问题方法+例题

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求边值问题例题

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程序)

(偏)微分方程一类边值问题的数值求解本文介绍了椭圆型(偏)微分方程一类边值问题的数值求解程序(笔者自编)及其使用方法。

程序基于差分原理,将连续的(偏)微分方程在节点上离散化,最终化为线性代数方程组。

对于原理方面的问题很多微分方程数值解的参考书中都有详尽的描述,但一般缺少具体的实现程序。

笔者认为,对这类数值方法是否理解并能运用的关键之一还是在于是否能写出用于解决问题的程序。

理论基础固然必不可少,程序确往往是我们解决问题的敲门砖。

一维椭圆型方程可表示如下:],[,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可以发现,程序仍能求解,但结果的光滑性不好。

偏微分方程边值问题的数值解法

偏微分方程边值问题的数值解法

求解偏微分方程的边值问题本实验学习使用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第一次上机作业

输入: >>tic, n=9;[u,k]=xsj(n), toc,surf(u)
计算Байду номын сангаас果如下:
u=
0 0 0 0 0 0 0 0 0 0 0 0 0.0256 0.0413 0.0508 0.0560 0.0577 0.0560 0.0508 0.0413 0.0256 0 0 0.0413 0.0686 0.0859 0.0955 0.0986 0.0955 0.0859 0.0686 0.0413 0 0.0508 0.0859 0.1088 0.1216 0.1258 0.1216 0.1088 0.0859 0.0508 0 0 0.0560 0.0955 0.1216 0.1364 0.1412 0.1364 0.1216 0.0955 0.0560 0 0 0.0577 0.0986 0.1258 0.1412 0.1462 0.1412 0.1258 0.0986 0.0577 0 0 0.0560 0.0955 0.1216 0.1364 0.1412 0.1364 0.1216 0.0955 0.0560 0 0 0.0508 0.0859 0.1088 0.1216 0.1258 0.1216 0.1088 0.0859 0.0508 0 0 0.0413 0.0686 0.0859 0.0955 0.0986 0.0955 0.0859 0.0686 0.0413 0 0 0.0256 0.0413 0.0508 0.0560 0.0577 0.0560 0.0508 0.0413 0.0256 0 0 0 0 0 0 0 0 0 0 0 0 0
程序一: function [u,k]=xsbj(n) % xsbj:用块 Jacobi 迭代法求解线性方程组 A*u=f % u:方程组的解; k 迭代次数; n:非边界点数 % a:方程组系数矩阵 Aii 的下对角线元素;b:方程系数矩阵 Aii 的主对角线元素 % c:方程组系数矩阵 Aii 的上对角线元素;d:追赶法所求方程的右端向量 % e:允许误差界;er:迭代误差 f=2*1/(n+1)^(2)*ones(n+2,n+2); a=-1*ones(1,n);b=4*ones(1,n);c=-1*ones(1,n);u=zeros(n+2,n+2);e=0.000000001; for k=1:2000 er=0; ub=u; for j=2:n+1 d(1:n)=f(2:n+1,j)+ub(2:n+1,j-1)+ub(2:n+1,j+1); x=zg(a,b,c,d); u(2:n+1,j)=x'; er=er+norm(ub(:,j)-u(:,j),1); end if er/n^2<e,break; end end 程序二: function x=zg(a,b,c,d) % zg:用追赶法求解线性方程组 n=length(b); % LU 分解 u(1)=b(1); for k=2:n if u(k-1)==0,D=0,return; end l(k)=a(k)/u(k-1); u(k)=b(k)-l(k)*c(k-1); end % 追赶法求解之追过程,求解 Ly=d y(1)=d(1); for k=2:n y(k)=d(k)-l(k)*y(k-1); end % 追赶法求解之追过程,求解 Ux=y if u(n)==0,D=0,return;end x(n)=y(n)/u(n); for k=n-1:-1:1 x(k)=(y(k)-c(k)*x(k+1))/u(k); end 输入:

matlab求解初边值问题的偏微分方程

matlab求解初边值问题的偏微分方程

偏微分方程是描述自然界中动态过程的重要数学工具,在工程领域中,求解偏微分方程是很多实际问题的重要一步。

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程序)

西京学数学软件实验任务书课程名称数学软件实验班级数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.有限差分法:基本思想是把连续的定解区域用有限个离散点构成的网格来代替,这些离散点称作网格的节点;把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似;把原方程和定解条件中的微商用差商来近似, 积分用积分和来近似,于是原微分方程和定解条件就近似地代之以代数方程组,即有限差分方程组 , 解此方程组就可以得到原问题在离散点上的近似解。

两点边值问题的不同迭代法比较及matlab实现

两点边值问题的不同迭代法比较及matlab实现

两点边值问题的不同迭代法比较及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求边值问题例题

matlab求边值问题例题
(原创实用版)
目录
1.MATLAB 求边值问题的基本概念
2.边值问题的例子
3.MATLAB 求解边值问题的步骤
4.总结
正文
一、MATLAB 求边值问题的基本概念
边值问题是数学物理中常见的问题,主要研究在边界上的函数值。

在MATLAB 中,我们可以通过构造方程组来求解边值问题。

边值问题可以应用于许多领域,如物理、工程和经济等。

二、边值问题的例子
假设有一个长方体,其中长、宽、高分别为 a、b、c,现在需要求解长方体内部任意一点的温度分布。

已知长方体的底部温度为 T1,顶部温度为 T2,左侧温度为 T3,右侧温度为 T4,底部和顶部的面积分别为 ab,长方体内部温度分布满足热传导方程。

这是一个典型的边值问题。

三、MATLAB 求解边值问题的步骤
1.构造方程组:根据问题的实际背景,构造相应的方程组。

例如,对于上述温度分布问题,我们可以构造以下方程组:
T/t = α(^2T/x^2 + ^2T/y^2 + ^2T/z^2)
T(x,y,0) = T1
T(x,y,c) = T2
T(0,y,z) = T3
T(a,y,z) = T4
2.设置边界条件:将边界条件加入方程组中。

3.求解方程组:使用 MATLAB 求解方程组,得到温度分布的数值解。

4.显示结果:使用 MATLAB 绘制温度分布的等值线图或三维图,直观地展示结果。

四、总结
MATLAB 作为一种强大的数学软件,可以方便地求解边值问题。

通过构造方程组、设置边界条件、求解方程组和显示结果等步骤,我们可以得到边值问题的数值解。

二阶椭圆偏微分方程实例求解(附matlab代码)

二阶椭圆偏微分方程实例求解(附matlab代码)

《微分方程数值解法》期中作业实验报告二阶椭圆偏微分方程第一边值问题姓名:学号:班级: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的微分方程边值问题数值解的算法研究_王来英

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带边值的微分方程求解

一、介绍MATLAB是一种强大的数学工具,能够用于解决各种数学问题,包括微分方程求解。

微分方程是描述动力系统、物理现象、工程问题等的重要数学工具,求解微分方程有助于我们了解和预测自然和人造系统的行为。

在实际问题中,经常会遇到带有边值条件的微分方程,对于这类问题,我们可以利用MATLAB进行求解。

二、MATLAB求解带边值的微分方程的基本步骤1. 定义微分方程我们需要将微分方程用MATLAB语言进行描述。

对于带边值的微分方程,我们通常会使用函数表示微分方程,并利用符号变量表示未知函数。

2. 离散化区域对于带边值的微分方程,通常需要在给定的区域上进行离散化处理。

这可以通过将区域划分为一系列离散点来实现。

在MATLAB中,我们可以使用网格生成函数或手动创建离散点进行离散化处理。

3. 设置边值条件在微分方程的求解过程中,需要考虑边值条件。

这些条件可以是方程在区域边界上的取值,也可以是方程在边界上的导数值。

在MATLAB中,我们可以使用条件语句或函数来设置边值条件。

4. 构建代数系统一旦我们定义了微分方程、离散化区域并设置了边值条件,就可以利用这些信息构建代数系统。

这可以通过将微分方程转化为离散点上的代数方程组来实现。

在MATLAB中,可以使用矩阵运算和线性方程求解函数来构建代数系统。

5. 求解代数系统我们可以使用MATLAB内置的求解器来求解代数系统,从而得到微分方程的近似解。

在求解之后,我们可以对得到的结果进行后处理,比如对解进行可视化或进行数值分析。

三、示例以一维波动方程为例,该方程可以表示为:∂^2u/∂t^2=c^2∂^2u/∂x^2在区间[0,1]上的边值条件为u(0,t)=0,u(1,t)=0,并且初始条件为u(x,0)=sin(πx),∂u/∂t(x,0)=0。

我们可以利用MATLAB对这个问题进行求解。

1. 定义微分方程我们可以使用MATLAB的符号计算工具箱定义波动方程,并引入符号变量。

matlab dsolve函数用法

matlab dsolve函数用法

matlab dsolve函数用法一、前言dsolve()函数是MATLAB中的一个常见工具,用于求解常微分方程初值问题和边值问题。

它可以方便地解决包括一阶和高阶常微分方程在内的的各类微分方程。

二、基本语法dsolve(func, var):其中,func是书写微分方程的字符串或匿名函数,变量var是代表未知函数的符号变量或字符串。

dsolve()函数可以自动识别微分方程的阶次,即使方程非线性、非均匀、变系数等复杂情形也没有问题。

同时,如果我们需要求解的是偏微分方程,也可以使用pdepe()函数。

三、简单示例下面给出一个简单的例子,首先定义一个一阶微分方程dy/dx=sin(x),然后求解其通解。

%定义微分方程syms y(x)ode = diff(y,x) == sin(x);%求解通解sol = dsolve(ode);%输出通解ySol(x) = sol;fprintf("通解为:\n");disp(ySol);代码解释:1. 首先定义一个符号变量y,它代表待求解的未知函数。

2. 然后定义这个函数的微分方程,这里使用syms定义符号变量,diff()函数表示取y(x)的导数。

3. 使用dsolve()函数解微分方程,并赋值给sol。

4. 最后将通解表示为一个函数ySol(x),并输出解。

输出结果:四、常见问题1. 存在多解如何求解?dsolve()函数默认输出一个解,但如果微分方程存在多个解的话,可以通过指定初始值以得到不同的解。

应该尽可能地将初始值选择在解的不同分支上,以获得不同的解。

2. 函数第二个参数可以是字符串吗?当我们输入变量名x时,我们必须在前面先使用syms声明符号变量x。

而如果我们直接输入一个字符串'x'作为函数的第二个参数,它是不需要使用syms声明的。

这就意味着,字符串类型的变量也是可以用作变量的。

3. dsolve()函数是否只适用于常微分方程?不是。

matlab 微分方程边界条件

matlab 微分方程边界条件

matlab 微分方程边界条件微分方程是数学中一个非常重要的概念,它描述了变量之间的关系及其变化规律。

在实际问题中,我们经常会遇到一些微分方程问题,而边界条件是解决微分方程问题中的重要一环。

本文将以边界条件为主题,探讨如何在MATLAB中处理微分方程问题,并给出一些实例进行说明。

我们来了解一下什么是边界条件。

在求解微分方程时,我们通常需要给出一些已知条件,这些条件被称为边界条件。

边界条件可以是函数在某些点处的值,也可以是函数在某些点处的导数值。

边界条件的作用是限定了微分方程的解的范围,使得求解问题具有唯一性。

在MATLAB中,求解微分方程可以使用ode45、ode23等数值方法。

这些方法可以解决一般的初值问题,但对于边值问题,即需要给定边界条件的情况,就需要使用其他方法。

MATLAB提供了一个专门用于求解边值问题的函数——bvp4c。

bvp4c函数的基本使用方法如下:```matlabsol = bvp4c(@(x,y)odefun(x,y),@(ya,yb)bcfun(ya,yb),guess);```其中,odefun是微分方程的函数句柄,bcfun是边界条件的函数句柄,guess是初值的猜测。

下面,我们通过一个具体的例子来说明如何使用bvp4c函数求解边值问题。

假设我们要求解如下的边值问题:$$y''(x) + y(x) = 0$$$$y(0) = 0, y(\pi) = 0$$我们需要将该二阶微分方程转化为一阶方程组的形式。

令$y_1(x) = y(x), y_2(x) = y'(x)$,则原方程可以转化为:$$y_1' = y_2$$$$y_2' = -y_1$$接下来,我们定义odefun函数和bcfun函数:```matlabfunction dydx = odefun(x,y)dydx = [y(2); -y(1)];endfunction res = bcfun(ya,yb)res = [ya(1); yb(1)];end```然后,我们可以调用bvp4c函数求解边值问题:```matlabsol = bvp4c(@(x,y)odefun(x,y),@(ya,yb)bcfun(ya,yb),guess);```其中,guess是初值猜测,可以根据实际情况进行调整。

matlab用差分法求解边值问题

matlab用差分法求解边值问题

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求边值问题例题

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. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

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,计算结果如何?
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];
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),
边值问题求解指令:bvp4c()
输出参数: 输出变量sol是一个结构体 sol.x是指令bvp4c所采用的网格节点; sol.y是y(x)在sol.x网点上的近似解值; sol.yp是y'(x)在sol.x网点上的近似解值; sol.parameters是微分方程所包含的未知参数的近似解值。 当被解微分方程包含未知参数时,该域存在。
z c z 0 c 1
z (0) 0, z (4) 2
令:
z1 z , z2 z1
边值问题的求解
求解:z ( 2q cos 2 x) z 0
q 15
边界条件:
solinit=bvpinit(linspace(0,pi,10),[1;1],lmb); opts=bvpset('Stats','on'); sol=bvp4c(@ODEfun,@BCfun,solinit,opts); lambda=sol.parameters x=[0:pi/60:pi]; y=deval(sol,x); plot(x,y(1,:),'b-',sol.x,sol.y(1,:),'ro') legend('解曲线','初始网格点解') % 定义ODEfun函数 function dydx=ODEfun(x,y,lmb) q=15; dydx=[y(2);-(lmb-2*q*cos(2*x))*y(1)]; % 定义BCfun函数 function bc=BCfun(ya,yb,lmb) bc=[ya(1)-1;ya(2);yb(2)];
边值问题求解指令: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等表示希望向被解微分方程传递的已知参数。如果无须 向微分方程传递参数,它们可以缺省。
使残差不断减小,从而获得满足精度要求的解
Matlab求解边值问题的基本指令
solinit=bvpinit(x,v,parameters) 生成bvp4c调用指令所必须的“解猜测网” sol=bvp4c(odefun,bcfun,solinit,options,p1,p2,…) 给出微分方程边值问题的近似解
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)];
z(0) 1, z (0) 0, z ( ) 0
本例中,微分方程与参数λ的数值有关。一般而言,对于任意的λ值,该问题无解, 但对于特殊的λ值(特征值),它存在一个解,这也称为微分方程的特征值问题。 对于此问题,可在bvpinit中提供参数的猜测值,然后重复求解BVP得到所需的 参数,返回参数为sol.parameters
边值问题的求解
y y 10 求解两点边值问题: y (0) y(1) 0
令:
y1 y , y2 y1
原方程组等价于以下标准形式的 方程组:
y y 1 2 y2 y1 10
边界条件为:
y1 (0) 0, y1 (1) 0
边值问题的求解
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];
相关文档
最新文档