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

李荣华二维抛物方程的初边值问题在数学和工程领域中扮演着重要角色。

通过使用Matlab,我们能够更深入地理解和解决这一问题。

在本文中,我们将从基础概念开始,逐步深入讨论李荣华二维抛物方程的初边值问题,并结合Matlab进行实际分析和解决。

一、初边值问题的概念和涵义1. 初边值问题的定义初边值问题是指在偏微分方程中,除了部分边界条件外,还需给出一些初始条件。

李荣华二维抛物方程的初边值问题即是在方程中给定初值条件和边界条件的问题。

2. 李荣华二维抛物方程简介李荣华二维抛物方程是描述热传导、扩散等现象的数学模型,在物理学和工程领域有着广泛的应用。

它的形式通常为一个关于未知函数和其在空间和时间上的导数的方程。

3. Matlab在初边值问题中的应用Matlab作为一种强大的数学建模和仿真工具,能够帮助我们求解各种类型的偏微分方程,包括初边值问题。

通过Matlab,我们可以对李荣华二维抛物方程进行数值求解和模拟,从而更清晰地理解问题本质。

二、李荣华二维抛物方程的初边值问题解析1. 方程的建立我们需要建立李荣华二维抛物方程的数学模型,并给定相应的初值条件和边界条件。

在Matlab中,我们可以通过定义矩阵和方程的形式来实现这一过程。

2. 数值求解我们利用Matlab提供的数值求解方法,如有限差分法或有限元法,对初边值问题进行求解。

通过程序的编写和运行,我们可以得到方程的数值解,并对物理过程有更直观的认识。

3. 结果分析在求得数值解之后,我们可以对结果进行分析和可视化。

通过Matlab 的绘图功能,我们能够直观地观察李荣华二维抛物方程的演化过程,以及初值条件和边界条件对解的影响。

三、个人观点和总结在本文中,我们深入探讨了李荣华二维抛物方程的初边值问题,并结合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有限元分析实例

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 微分方程边界条件

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有限元分析实例
唯一的一个函数stima 用来计算各单元刚度矩阵。为考虑进一步扩展之可能性,我们采用了等参单元的表达形式,引入了雅可比矩阵 J。事实上,该函数可直接用于三维空间四面体单元刚度矩阵之计算。
求解器之使用
可以使用文本文件存储输入数据,用MATLAB 提供的 load 命令载入。仅需下面三行代码即可完成读取输入数据,调用求解器求解,绘制结果云图诸个步骤。
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);
输入数据格式:
nds:节点列表矩阵。对于 [公式] 个节点,为一 [公式] 矩阵。第一列为 [公式] 坐标值,第二列为 [公式] 坐标值。
其中[公式] 为检验函数,[公式] 为一适当索伯列夫空间(在本例中也是一希尔伯特空间)。[公式] 为一双线性型,[公式] 为一线性型。其具体表达式为
[公式]
有限元空间离散
我们采用最简单的二维[公式] 单元离散 [公式]。[公式] 单元即三节点线性三角形单元,其插值基函数(即形函数) [公式] 为一次多项式。有限元之核心思想为:使用离散的函数空间 [公式] 来分片逼近连续的函数空间 [公式]。于是所求之近似解 [公式] 可写成基函数之线性组合
[公式]
其中[公式] 为待求系数,常称为自由度。将此近似解之表达式代入前述问题之弱形式,并取检验函数 [公式],可得
[公式]
写成矩阵形式,便成为常见之有限元方程
[公式]
由于历史原因,通常采用固体力学中的习惯命名:[公式] 为刚度矩阵,[公式] 为自由度向量,[公式] 为载荷向量。

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)。

初始解生成函数:bvpinit() 初始解生成函数 p ()
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求解边值问题的基本指令
solinit=bvpinit(x,v,parameters) 生成bvp4c b 4 调用指令所必须的“解猜测网” sol=bvp4c(odefun,bcfun,solinit,options,p1,p2,…) sol=bvp4c(odefun bcfun solinit options p1 p2 ) 给出微分方程边值问题的近似解 sxint=deval(sol,xint) 计算微分方程积分区间内任何一点的解值
z (0) 1, z (0) 0, z ( ) 0
本例中,微分方程与参数λ的数值有关。一般而言,对于任意的λ值,该问题无解, 但对于特殊的λ值(特征值),它存在一个解,这也称为微分方程的特征值问题。 对于此问题,可在bvpinit中提供参数的猜测值,然后重复求解BVP得到所需的 参数,返回参数为sol.parameters
边值问题的求解
如果取1,计算结果如何?
infinity=6; solinit=bvpinit(linspace(0 infinity 5) [0 0 1]); solinit=bvpinit(linspace(0,infinity,5),[0 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 fprintf( Value computed here is f''''(0)=%7 f (0)=%7.5f.\n 5f \n',f(3,1)) 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];
边值问题的求解
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]; x=[0:0.1:1]; y=deval(sol,x); xP=[0:0.1:1.0]; yP=[1 y [ 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') Northwest ) 'location' 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]; bc=[ya(1)-1;yb(1)-3];
边值问题求解指令:bvp4c() 边值问题求解指令 p ()
输出参数: 输出变量sol是一个结构体 sol.x是指令bvp4c所采用的网格节点; sol.y是y(x)在sol.x网点上的近似解值; sol.yp sol yp是y'(x)在sol.x sol x网点上的近似解值; 网点上的近似解值 sol.parameters是微分方程所包含的未知参数的近似解值。 当被解微分方程包含未知参数时,该域存在。
边值问题求解指令:bvp4c() 边值问题求解指令 p ()
sol=bvp4c(odefun,bcfun,solinit,options,p1,p2,…) p ( , , , p ,p ,p , ) 输入参数:
odefun是计算导数的m函数文件。该函数的基本形式为: dydx=odefun(x y parameters p1 p2 ),在此,自变量 dydx=odefun(x,y,parameters,p1,p2,…) 在此 自变量x是标量, 是标量 y, dydx是列向量。 bcfun是计算边界条件下残数的m函数文件。其基本形式为: res=bcfun(ya,yb,parameters,p1,p2,…) b f ( b t 1 2 ),文件输入宗量 文件输入宗量ya,yb b是边界 条件列向量,分别代表y在a和b处的值。res是边界条件满足时的残数 列向量。注意:例如odefun函数的输入宗量中包含若干“未知”和 “已知”参数,那么不管在边界条件计算中是否用到,它们都应作为 bcfun的输入宗量。 p 是用来改变bvp4c p 算法的控制参数的。在最基本用法 输入宗量options 中,它可以缺省,此时一般可以获得比较满意的边值问题解。如需更 改可采用bvpset函数,使用方法同odeset函数。 输入宗量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); 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=ODEfun(x y) dydx=[y(2);y(1)+10]; % 定义BCfun函数 function bc=BCfun(ya,yb) bc=[ya(1);yb(1)];
边值问题的求解
y y 10 求解两点边值问题令:
y1 y, y2 y1
原方程组等价于以下标准形式的 方程组:
y1 y2 y2 y1 10
边界条件为:
y1 (0) 0, 0 y1 (1) 0
y1 y, y2 y1
原方程组等价于以下标准形式的 方程组:
y1 y2 2 y x (1 ) y1 1 2
边界条件为:
y1 (0) 1 0 y1 (1) 3 0
边值问题的求解
求解:
c=1; solinit=bvpinit(linspace(0,4,10),[1 1]); sol=bvp4c(@ODEfun sol bvp4c(@ODEfun,@BCfun,solinit,[],c); @BCfun solinit [] c); x=[0:0.1:4]; y=deval(sol,x); plot(x,y(1,:),'b-',sol.x,sol.y(1,:),'ro') plot(x,y(1,:), b ,sol.x,sol.y(1,:), ro ) legend('解曲线','初始网格点解') % 定义ODEfun函数 function dydx=ODEfun(x,y,c) y ( ,y, ) dydx=[y(2);-c*abs(y(1))]; % 定义BCfun函数 function bc=BCfun(ya,yb,c) bc=[ya(1);yb(1)+2];
工程应用中我们经常遇到边值问题,这些是那些ode**函数无能为力的, 当然我们可以自己编写函数求解(shooting),但是不是能力所及的. Matlab中提供了 bvp解算器。 解算器 solinit = bvpinit(x, yinit, params) sol = bvpsolver(odefun,bcfun,solinit,options) 由于边值问题可能有多解,为了便于我们确定那个解是我们需要的,所 以必须使用bvpinit函数对初值进行估计 解算器(bvpsolver) (b l ): Matlab中提供了bvp4c和bvp5c,后者误差控制更好些 输 参数 输入参数: x:需要计算的网格点,相当于ode**的tspan yinit:猜测的值,可以是具体值,也可以是函数,类似与 ode**的 x0 params:其它未知参数,也是一个猜测值 其它未知参数 也是 个猜测值 odefun:描述边值问题微分方程的函数句柄 (x 的上下限即认为两个边界),支持多 bcfun:边值函数,一般是双边值( 边值 solinit:由bvpinit生成的初始化网格 options:BVP解算器优化参数,可以通过 解算器优化参数 可以通过bvpset设置 输出参数:大部分同理ode45
相关文档
最新文档