偏微分方程数值解法MATLAB源码
matlab求解最简单的一阶偏微分方程
matlab求解最简单的一阶偏微分方程一、引言在科学和工程领域，偏微分方程是非常重要的数学工具，用于描述各种现象和过程。
而MATLAB作为一种强大的数值计算软件，可以用来求解各种复杂的偏微分方程。
本文将以MATLAB求解最简单的一阶偏微分方程为主题，探讨其基本原理、数值求解方法以及具体实现过程。
二、一阶偏微分方程的基本原理一阶偏微分方程是指只含有一个未知函数的偏导数的微分方程。
最简单的一阶偏微分方程可以写成如下形式：\[ \frac{\partial u}{\partial t} = F(x, t, u, \frac{\partial u}{\partial x}) \]其中，\(u(x, t)\) 是未知函数，\(F(x, t, u, \frac{\partial u}{\partial x})\) 是给定的函数。
一阶偏微分方程可以描述很多实际问题，比如热传导、扩散等。
在MATLAB中，我们可以使用数值方法求解这类方程。
三、数值求解方法1. 有限差分法有限差分法是一种常用的数值求解偏微分方程的方法。
其基本思想是用离散的方式来逼近偏导数，然后将偏微分方程转化为代数方程组。
在MATLAB中，我们可以使用内置的求解器来求解离散化后的代数方程组。
2. 特征线法特征线法是另一种常用的数值求解方法，它利用特征线方程的特点来求解偏微分方程。
这种方法在求解一维情况下的偏微分方程时特别有效，可以提高求解的效率和精度。
四、MATLAB求解过程在MATLAB中，我们可以使用`pdepe`函数来求解一阶偏微分方程。
该函数可以针对特定的方程和边界条件，利用有限差分法进行离散化求解。
下面给出一个具体的例子来说明如何使用MATLAB求解最简单的一阶偏微分方程。
假设我们要求解如下的一维热传导方程：\[ \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} \]其中，\(\alpha\) 是热传导系数。
Matlab求解微分方程(组)及偏微分方程(组)
第四道 Matlab供解微分圆程（组）之阳早格格创做表里介绍：Matlab供解微分圆程（组）下令供解真例：Matlab供解微分圆程（组）真例本质应用问题通过数教修模所归纳得到的圆程，绝大普遍皆是微分圆程，真真能得到代数圆程的机会很少.另一圆里，不妨供解的微分圆程也是格中有限的，特地是下阶圆程战偏偏微分圆程（组）.那便央供咱们必须钻研微分圆程（组）的解法：剖析解法战数值解法.一．相关函数、下令及简介1.正在Matlab中，用大写字母D表示导数，Dy表示y关于自变量的一阶导数，D2y表示y关于自变量的二阶导数，依此类推.函数dsolve用去办理常微分圆程（组）的供解问题，调用要领为：X=dsolve(‘eqn1’,’eqn2’,…)函数dsolve用去解标记常微分圆程、圆程组，如果不初初条件，则供出通解，如果有初初条件，则供出特解.注意，系统缺省的自变量为t2.函数dsolve供解的是常微分圆程的透彻解法，也称为常微分圆程的标记解.然而是，有洪量的常微分圆程虽然从表里上道，其解是存留的，然而咱们却无法供出其剖析解，此时，咱们需要觅供圆程的数值解，正在供常微分圆程数值解圆里，MATLAB 具备歉富的函数，咱们将其统称为solver ，其普遍要领为：[T,Y]=solver(odefun,tspan,y0)证明：(1)solver 为下令ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 、ode15i 之一.(2)odefun 是隐现微分圆程'(,)y f t y =正在积分区间tspan 0[,]f t t =上从0t 到f t 用初初条件0y 供解.(3)如果要赢得微分圆程问题正在其余指定时间面012,,,,f t t t t 上的解，则令tspan 012[,,,]f t t t t =(央供是单调的).(4)果为不一种算法不妨灵验的办理所有的ODE 问题，为此，Matlab 提供了多种供解器solver ，对付于分歧的ODE 问题，采与分歧的solver.表1 Matlab 华文本文献读写函数证明：ode23、ode45是极其时常使用的用去供解非刚刚性的尺度形式的一阶微分圆程（组）的初值问题的解的Matlab时常使用步调，其中：ode23采与龙格-库塔2阶算法，用3阶公式做缺面预计去安排步少，具备矮等的粗度.ode45则采与龙格-库塔4阶算法，用5阶公式做缺面预计去安排步少，具备中等的粗度.3．正在matlab下令窗心、步调或者函数中创修局部函数时，可用内联函数inline，inline函数形式相称于编写M 函数文献，然而不需编写M-文献便不妨形貌出某种数教关系.调用inline函数，只可由一个matlab表白式组成，而且只可返回一个变量，不允许[u,v]那种背量形式.果而，所有央供逻辑运算或者乘法运算以供得最后截止的场合，皆不克不迭应用inline函数，inline函数的普遍形式为：FunctionName=inline(‘函数真质’, ‘所有自变量列表’)比圆：（供解F(x)=x^2*cos(a*x)-b ,a,b是标量；x是背量）正在下令窗心输进:Fofx=inline(‘x .^2*cos(a*x)-b’ , ‘x’,’a’,’b’);g= Fofx([pi/3 pi/3.5],4,1)注意：由于使用内联对付象函数inline 不需要其余修坐m 文献，所有使用比较便当，其余正在使用ode45函数的时间，定义函数往往需要编写一个m 文献去单独定义，那样便当于管造文献，那里不妨使用inline 去定义函数.二．真例介绍例1 供解微分圆程2'2x y xy xe -+= 步调：syms x y; y=dsolve(‘Dy+2*x*y=x*exp(-x^2)’,’x ’) 例2 供微分圆程'0x xy y e +-=正在初初条件(1)2y e =下的特解并画出解函数的图形.步调：syms x y; y=dsolve(‘x*Dy+y-exp(1)=0’,’y(1)=2*exp(1)’,’x ’);ezplot(y)例 3 供解微分圆程组530t dx x y e dt dy x y dt ⎧++=⎪⎪⎨⎪--=⎪⎩正在初初条件00|1,|0t t x y ====下的特解并画出解函数的图形.步调：syms x y t[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t')simple(x);simple(y)ezplot(x,y,[0,1.3]);axis auto2.用ode23、ode45等供解非刚刚性尺度形式的一阶微分圆程（组）的初值问题的数值解（近似解）例4 供解微分圆程初值问题2222(0)1dy y x x dx y ⎧=-++⎪⎨⎪=⎩的数值解，供解范畴为区间[0,0.5].步调：fun=inline('-2*y+2*x^2+2*x','x','y');[x,y]=ode23(fun,[0,0.5],1);plot(x,y,'o-')例5 供解微分圆程22'2(1)0,(0)1,(0)0d y dy y y y y dt dt μ--+===的解，并画出解的图形.12,,7dy x y x dt μ===，则编写M-文献function fy=vdp(t,x)fy=[x(2);7*(1-x(1)^2)*x(2)-x(1)];end正在Matlab 下令窗心编写步调y0=[1;0][t,x]=ode45(@vdp,[0,40],y0);或者[t,x]=ode45('vdp',[0,40],y0);y=x(:,1);dy=x(:,2);plot(t,y,t,dy)训练与思索：M-文献vdp.m 改写成inline 函数步调？Euler 合线法供解的基础思维是将微分圆程初值问题 化成一个代数(好分)圆程，主要步调是用好商()()y x h y x h +-代替微商dydx ，于是 记1,(),k k k k x x h y y x +=+=进而1(),k k y y x h +=+于是例6用Euler 合线法供解微分圆程初值问题的数值解（步少h 与），供解范畴为区间[0,2].分解：本问题的好分圆程为步调：>> clear>> f=sym('y+2*x/y^2');>> a=0;>> b=2;>> h=0.4;>> n=(b-a)/h+1;>> x=0;>> y=1;>> szj=[x,y];%数值解>> for i=1:n-1y=y+h*subs(f,{'x','y'},{x,y});%subs ，替换函数x=x+h;szj=[szj;x,y];end>>szj>> plot(szj(:,1),szj(:,2))证明：替换函数subs比圆：输进subs(a+b,a,4) 意义便是把a用4替换掉，返回4+b，也不妨替换多个变量，比圆：subs(cos(a)+sin(b),{a,b},[sym('alpha'),2])分别用字符alpha替换a战2替换b，返回 cos(alpha)+sin(2)特地证明：本问题可进一步利用四阶Runge-Kutta法供解，Euler合线法本质上便是一阶Runge-Kutta法，Runge-Kutta法的迭代公式为相映的Matlab步调为：>> clear>> f=sym('y+2*x/y^2');>> a=0;>> b=2;>> h=0.4;>> n=(b-a)/h+1;>> x=0;>> y=1;>> szj=[x,y];%数值解>> for i=1:n-1l1=subs(f,{'x','y'},{x,y});替换函数l2=subs(f,{'x','y'},{x+h/2,y+l1*h/2});l3=subs(f,{'x','y'},{x+h/2,y+l2*h/2});l4=subs(f,{'x','y'},{x+h,y+l3*h});y=y+h*(l1+2*l2+2*l3+l4)/6;x=x+h;szj=[szj;x,y];end>>szj>> plot(szj(:,1),szj(:,2))训练与思索：(1)ode45供解问题并比较好别.(2)利用Matlab 供微分圆程(4)(3)''20y y y -+=的解. (3)供解微分圆程''2',2(1)0,030,(0)1,(0)0y y y y x y y --+=≤≤==的特解.(4)利用Matlab 供微分圆程初值问题2''''00(1)2,|1,|3x x x y xy y y ==+===的解.指示：尽大概多的思量解法三．微分圆程变更为一阶隐式微分圆程组Matlab 微分圆程解算器只可供解尺度形式的一阶隐式微分圆程（组）问题，果此正在使用ODE 解算器之前，咱们需要干的第一步，也是最要害的一步便是借帮状态变量将微分圆程（组）化成Matlab 可交受的尺度形式.天然，如果ODEs 由一个或者多个下阶微分圆程给出，则咱们应先将它变更成一阶隐式常微分圆程组.底下咱们以二个下阶微分圆程组形成的ODEs 为例介绍怎么样将它变更成一个一阶隐式微分圆程组.Step 1 将微分圆程的最下阶变量移到等式左边，其余移到左边，并按阶次从矮到下排列.形式为：Step 2 为每一阶微分式采用状态变量，最下阶除中注意：ODEs 中所有是果变量的最下阶次之战便是需要的状态变量的个数，最下阶的微分式不需要给它状态变量.Step 3 根据采用的状态变量，写出所有状态变量的一阶微分表白式训练与思索：(1)供解微分圆程组 其中2r =1r =*1,μμ=-1/82.45,μ=(0) 1.2,x =(2)供解隐式微分圆程组提示：使用标记预计函数solve供'''',x y ，而后利用供解微分圆程的要领四．偏偏微分圆程解法Matlab 提供了二种要领办理PDE 问题，一是使用pdepe 函数，它不妨供解普遍的PDEs,具备较大的通用性，然而只收援下令形式调用；二是使用PDE 工具箱，不妨供解特殊PDE 问题，PDEtoll 有较大的限造性，比圆只可供解二阶PDE 问题，而且不克不迭办理片微分圆程组，然而是它提供了GUI 界里，从搀纯的编程中解脱出去，共时还不妨通过File —>Save As 直交死成M 代码.1.普遍偏偏微分圆程（组）的供解(1)Matlab 提供的pdepe 函数，不妨直交供解普遍偏偏微分圆程（组），它的调用要领为：sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)@pdefun 是PDE 的问题形貌函数，它必须换成尺度形式：那样，PDE 便不妨编写出心函数：[c,f,s]=pdefun(x,t,u,du)，m,x,t 对付应于式中相关参数，du 是u 的一阶导数，由给定的输进变量可表示出c,f,s 那三个函数.@pdebc 是PDE 的鸿沟条件形貌函数，它必须化为形式：于是边值条件不妨编写函数形貌为：[pa,qa,pb,qb]=pdebc(x,t,u,du),其中a 表示下鸿沟，b 表示上鸿沟.@pdeic 是PDE 的初值条件，必须化为形式：00(,)u x t u =，故不妨使用函数形貌为：u0=pdeic(x)sol 是一个三维数组，sol(:,:,i)表示i u 的解，换句话道，k u 对付应x(i)战t(j)时的解为sol(i,j,k)，通过sol ，咱们不妨使用pdeval 函数直交预计某个面的函数值.(2)真例证明供解偏偏微分其中， 5.7311.46()x x F x e e -=-且谦脚初初条件12(,0)1,(,0)0u x u x ==及鸿沟条件1(0,)0,u t x ∂=∂221(0,)0,(1,)1,(1,)0u u t u t t x ∂===∂解：(1)对付照给出的偏偏微分圆程战pdepe 函数供解的尺度形式，本圆程改写为 可睹1121220.024()10,,,()10.17u F u u x m c f s F u u u x ∂⎡⎤⎢⎥--⎡⎤⎡⎤∂====⎢⎥⎢⎥⎢⎥-∂⎣⎦⎣⎦⎢⎥⎢⎥∂⎣⎦%目标PDE 函数function [c,f,s]=pdefun(x,t,u,du)c=[1;1];f=[0.024*du(1);0.17*du(2)];temp=u(1)-u(2);s=[-1;1].*(exp(5.73*temp)-exp(-11.46*temp))end(2)鸿沟条件改写为：下鸿沟2010.*00f u ⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦上鸿沟1110.*000u f -⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ %鸿沟条件函数function [pa,qa,pb,qb]=pdebc(xa,ua,xb,ub,t)pa=[0;ua(2)];qa=[1;0];pb=[ub(1)-1;0];qb=[0;1];end(3)初值条件改写为：1210u u ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦%初值条件函数function u0=pdeic(x)u0=[1;0];end(4)编写主调函数clcx=0:0.05:1;t=0:0.05:2;m=0;sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t);subplot(2,1,1)surf(x,t,sol(:,:,1))subplot(2,1,2)surf(x,t,sol(:,:,2))训练与思索： This example illustrates the straightforward formulation, computation, and plotting of the solution of a single PDE.This equation holds on an interval 01x ≤≤ for times 0t ≥. The PDE satisfies the initial condition (,0)sin u x x π= and boundary conditions(1)PDEtool （GUI ）供解偏偏微分圆程的普遍步调正在Matlab 下令窗心输进pdetool ，回车，PDE 工具箱的图形用户界里(GUI)系统便开用了.从定义一个偏偏微分圆程问题到完毕解偏偏微分圆程的定解，所有历程大概不妨分为六个阶段Step 1 “Draw模式”画造仄里有界地区Ω，通过公式把Matlab系统提供的真体模型：矩形、圆、椭圆战多边形，拉拢起去，死成需要的仄里地区.Step 2 “Boundary模式”定义鸿沟，声明分歧鸿沟段的鸿沟条件.Step 3 “PDE模式”定义偏偏微分圆程，决定圆程典型战圆程系数c,a,f,d，根据简直情况，还不妨正在分歧子地区声明分歧系数.Step 4 “Mesh模式”网格化地区Ω，不妨统造自动死成网格的参数，对付死成的网格举止多次细化，使网格分隔更细更合理.Step 5 “Solve模式”解偏偏微分圆程，对付于椭圆型圆程不妨激活并统造非线性自符合解题器去处理非线性圆程；对付于扔物线型圆程战单直型圆程，树坐初初鸿沟条件后不妨供出给定时刻t的解；对付于特性值问题，不妨供出给定区间上的特性值.供解完毕后，不妨返回到Step 4，对付网格进一步细化，举止再次供解.Step 6 “View模式”预计截止的可视化，不妨通过树坐系统提供的对付话框，隐现所供的解的表面图、网格图、等下线图战箭头梯形图.对付于扔物线型战单直线型问题的解还不妨举止径画演示.(2)真例证明用法供解一个正圆形地区上的特性值问题：正圆形地区为：11,1 1.-≤≤-≤≤x x(1)使用PDE工具箱挨开GUI供解圆程(2)加进Draw模式，画造一个矩形，而后单打矩形，正在弹出的对付话框中树坐Left=-1,Bottom=-1,Width=2,Height=2,确认并关关对付话框(3)加进Boundary模式，鸿沟条件采与Dirichlet条件的默认值(4)加进PDE模式，单打工具栏PDE按钮，正在弹出的对付话框中圆程典型采用Eigenmodes,参数树坐c=1,a=-1/2,d=1,确认后关关对付话框(5)单打工具栏的∆按钮，对付正圆形地区举止初初网格剖分，而后再对付网格进一步细化剖分一次(6)面开solve菜单，单打Parameters选项，正在弹出的对付话框中树坐特性值地区为[-20,20](7)单打Plot菜单的Parameters项，正在弹出的对付话框中选中Color、Height(3-D plot)战show mesh项，而后单打Done确认(8)单打工具栏的“=”按钮，开初供解。
(偏)微分方程一类边值问题的数值求解(附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实验报告
偏微分方程数值实验报告八实验题目：利用有限差分法求解.0)1(,0)1(),()()(==-=+''-u u x f x u x u 真解为)1()(22x e x u x -=-实现算法：对于两点边值问题,)(,)(,,d 22βα==∈=-b u a u l x f dxu(1)其中),(b a l =f b a ),(<为],[b a l =上的连续函数，βα，为给定常数.其相应的有限差分法的算法如下：1.对求解区域做网格剖分，得到计算网格.在这里我们对区间l 均匀剖分n 段，每个剖分单元的剖分步长记为nab h -=.2.对微分方程中的各阶导数进行差分离散，得到差分方程.运用的离散方法有：方法一:用待定系数和泰勒展开进行离散)()()()(d )(d 111122++--++≈i i i i i i i i x u x u x u x x u ααα方法二：利用差商逼近导数21122)()(2)()(d )(d h x u x u x u x x u i i i i i -++-≈(2)将(2)带入(1)可以得到)()()()(2)(211u R x f h x u x u x u i i i i i +=+---+,其中)(u R i 为无穷小量，这时我们丢弃)(u R i ，则有在i x 处满足的计算公式：1,...,1)()()(2)(211-==+---+n i x f hx u x u x u i i i i ，(3)3.根据边界条件，进行边界处理.由(1)可得βα==n u u ,0(4)称(3)(4)为逼近(1)的差分方程，并称相应的数值解向量1-n U 为差分解，i u 为)(i x u 的近似值.4.最后求解线性代数方程组，得到数值解向量1-n U .实验题目：用Lax-Wendroff 格式求解方程：.4sin 1),1(],1,0[,2sin 1)0,(,0),1,0(,02t t u x x x u t x u u x t ππ+=∈+=>∈=- （1） （精确解).2(2sin 1t x u ++=π） 数值边值条件分别为： （a ）).(20101n 0nn nu u hu u -+=+τ （b ）.1n 0nu u =（c ）.02-12111n 0=++++n n u u u请将计算结果与精确解进行比较。
matlab解偏微分方程
ui,j +1 − ui,j = Hui,j +1 ∆t Hui,j = a2 ui+1,j − 2ui,j + ui−1,j (∆x)2
ui,j +1 − ui,j = Hui,j ∆t 将显式与隐式相加,得平均公式 ui,j +1 − ui,j 1 1 = Hui,j + Hui,j +1 ∆t 2 2
得ui,0 = ui,2 − 2ψi
1 ui,2 = [c(ui+1,1 + ui−1,1) + 2(1 − c)ui,1 + 2ψi t] 2
3.3
例题 两端固定的弦振动
两端固定的弦, 初速为零，初位移是 h x, (0 ≤ x ≤ 2/3) 2 / 3 u(x, 0) = 1−x , (2/3 < x ≤ 1) h 1 − 2/3
作图所用程序如下,其中取c = 0.05, l = 1, h = 0.05.这里使用的方程 与初始条件表示方法与上一节相同. N=4000; c=0.05; x=linspace(0,1,420)’; u1(1:420)=0; u2(1:420)=0; u3(1:420)=0; u1(2:280)=0.05/279*(1:279)’; u1(281:419)=0.05/(419-281)*(419-(281:419)’); u2(2:419)=u1(2:419)+c/2*(u1(3:420)-2*u1(2:419)+u1(1:418)); h=plot(x,u1,’linewidth’,3); axis([0,1,-0.05,0.05]); set(h,’EraseMode’,’xor’,’MarkerSize’,18) for k=2:N set(h,’XData’,x,’YData’,u2) ; drawnow; u3(2:419)=2*u2(2:419)-u1(2:419)+c*(u2(3:420)... -2*u2(2:419)+u2(1:418)); u1=u2; u2=u3; end
matlab、lingo程序代码21-微分方程【精选】
MATLAB 求微分方程的解一、问题背景与实验目的实际应用问题通过数学建模所归纳而得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少．另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）．这就要求我们必须研究微分方程（组）的解法，既要研究微分方程（组）的解析解法（精确解），更要研究微分方程（组）的数值解法（近似解）．对微分方程（组）的解析解法(精确解)，Matlab 有专门的函数可以用，本实验将作一定的介绍．本实验将主要研究微分方程(组)的数值解法（近似解），重点介绍 Euler 折线法．二、相关函数（命令）及简介1．dsolve('equ1','equ2',…)：Matlab 求微分方程的解析解．equ1、equ2、…为方程（或条件）．写方程（或条件）时用 Dy 表示y 关于自变量的一阶导数，用用 D2y 表示 y 关于自变量的二阶导数，依此类推．2．simplify(s)：对表达式 s 使用 maple 的化简规则进行化简．例如：syms xsimplify(sin(x)^2 + cos(x)^2)ans=13．[r,how]=simple(s)：由于 Matlab 提供了多种化简规则，simple 命令就是对表达式 s 用各种规则进行化简，然后用 r 返回最简形式，how 返回形成这种形式所用的规则．例如：syms x[r,how]=simple(cos(x)^2-sin(x)^2)r = cos(2*x)how = combine4．[T,Y] = solver(odefun,tspan,y 0) 求微分方程的数值解．说明：(1) 其中的 solver 为命令 ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 之一．(2) odefun 是显式常微分方程：⎪⎩⎪⎨⎧==00)(),(y t y y t f dtdy(3) 在积分区间 tspan =上，从到，用初始条件求解．],[0f t t 0t f t 0y (4)要获得问题在其他指定时间点上的解，则令tspan =,210,,t t t（要求是单调的）．],,,[,210f t t t t (5) 因为没有一种算法可以有效地解决所有的 ODE 问题，为此，Matlab 提供了多种求解器 Solver ，对于不同的ODE 问题，采用不同的Solver ．求解器Solver ODE 类型特点说明ode45非刚性单步算法；4、5阶Runge-Kutta 方程；累计截断误差达3)(x ∆大部分场合的首选算法ode23非刚性单步算法；2、3阶Runge-Kutta 方程；累计截断误差达3)(x ∆使用于精度较低的情形ode113非刚性多步法；Adams 算法；高低精度均可到6310~10--计算时间比 ode45 短ode23t 适度刚性采用梯形算法适度刚性情形ode15s刚性多步法；Gear's 反向数值微分；精度中等若 ode45 失效时，可尝试使用ode23s刚性单步法；2阶 Rosebrock 算法；低精度当精度较低时，计算时间比 ode15s 短ode23tb刚性梯形算法；低精度当精度较低时，计算时间比 ode15s 短(6) 要特别的是：ode23、ode45 是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的 Matlab 的常用程序，其中：ode23 采用龙格-库塔2 阶算法，用3 阶公式作误差估计来调节步长，具有低等的精度．ode45 则采用龙格-库塔4 阶算法，用5 阶公式作误差估计来调节步长，具有中等的精度．5．ezplot(x,y,[tmin,tmax])：符号函数的作图命令．x,y 为关于参数t 的符号函数，[tmin,tmax] 为 t 的取值范围．6．inline()：建立一个内联函数．格式：inline('expr', 'var1', 'var2',…) ，注意括号里的表达式要加引号．例：Q = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)三、实验内容1. 几个可以直接用 Matlab 求微分方程精确解的例子：例1：求解微分方程，并加以验证．22x xe xy dxdy-=+求解本问题的Matlab 程序为：1syms x y %line1y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') %line2diff(y,x)+2*x*y-x*exp(-x^2) %line3simplify(diff(y,x)+2*x*y-x*exp(-x^2)) %line4说明：(1) 行line1是用命令定义x,y 为符号变量．这里可以不写，但为确保正确性，建议写上；(2) 行line2是用命令求出的微分方程的解：1/2*exp(-x^2)*x^2+exp(-x^2)*C1(3) 行line3使用所求得的解．这里是将解代入原微分方程，结果应该为0，但这里给出：-x^3*exp(-x^2)-2*x*exp(-x^2)*C1+2*x*(1/2*exp(-x^2)*x^2+exp(-x^2)*C1)(4) 行line4 用 simplify() 函数对上式进行化简，结果为 0， 表明的确是微分方程的解．)(x y y =例2：求微分方程在初始条件下的特解，并画出解函0'=-+x e y xy e y 2)1(=数的图形．求解本问题的 Matlab 程序为：syms x yy=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x')ezplot(y)微分方程的特解为：y=1/x*exp(x)+1/x* exp (1) (Matlab 格式)，即，xe e y x+=解函数的图形如图 1：图1例3：求微分方程组在初始条件下的特解，⎪⎪⎩⎪⎪⎨⎧=--=++035y x dtdy e y x dtdx t0|,1|00====t t y x 课，系列话并画出解函数的图形．求解本问题的 Matlab 程序为：syms x y t[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t')simple(x);simple(y);ezplot(x,y,[0,1.3]);axis auto微分方程的特解(式子特别长)以及解函数的图形均略．2. 用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)．例4：求解微分方程初值问题的数值解，求解范围为⎪⎩⎪⎨⎧=++-=1)0(2222y xx y dx dy 区间[0, 0.5]．fun=inline('-2*y+2*x^2+2*x','x','y');[x,y]=ode23(fun,[0,0.5],1);x';y';plot(x,y,'o-')>> x'ans =0.0000 0.0400 0.0900 0.1400 0.1900 0.24000.2900 0.3400 0.3900 0.4400 0.4900 0.5000>> y'ans =1.0000 0.9247 0.8434 0.7754 0.7199 0.67640.6440 0.6222 0.6105 0.6084 0.6154 0.6179图形结果为图 2．图2例 5：求解描述振荡器的经典的 Ver der Pol 微分方程.7,0)0(',1)0(,0)1(222====+--μμy y y dt dy y dt y d 分析：令则,,121dt dx x y x ==.)1(,1221221x x x dtdx x dt dx --==μ先编写函数文件verderpol.m ：function xprime = verderpol(t,x)global mu;xprime = [x(2);mu*(1-x(1)^2)*x(2)-x(1)];再编写命令文件vdp1.m ：global mu;mu = 7;y0=[1;0][t,x] = ode45('verderpol',[0,40],y0);x1=x(:,1);x2=x(:,2);plot(t,x1)图形结果为图3．图33. 用 Euler 折线法求解前面讲到过，能够求解的微分方程也是十分有限的．下面介绍用 Euler 折线法求微分方程的数值解(近似解)的方法．Euler 折线法求解的基本思想是将微分方程初值问题⎪⎩⎪⎨⎧==00)(),,(y x y y x f dx dy化成一个代数方程，即差分方程，主要步骤是用差商替代微商hx y h x y )()(-+ 党，做习贯彻严三实创先,于是：dxdy⎪⎩⎪⎨⎧==-+)()),(,()()(00x y y x y x f h x y h x y k k k k 记，从而，则有)(,1k k k k x y y h x x =+=+)(1h x y y k k +=+1,,2,1,0).,(,),(1100-=⎪⎩⎪⎨⎧+=+==++n k y x hf y y h x x x y y k k k k k k 例 6：用 Euler 折线法求解微分方程初值问题⎪⎩⎪⎨⎧=+=1)0(,22y y x y dxdy 的数值解(步长h 取0.4)，求解范围为区间[0,2]．解：本问题的差分方程为1,,2,1,0).2),( ),(,,4.0,1,021100-=⎪⎪⎪⎩⎪⎪⎪⎨⎧+=+=+====++n k y x y y x f y x hf y y h x x h y x k k k k k k （其中：相应的Matlab 程序见附录 1．数据结果为： 0 1.0000 0.4000 1.4000 0.8000 2.1233 1.2000 3.1145 1.6000 4.4593 2.0000 6.3074图形结果见图4：图4特别说明：本问题可进一步利用四阶 Runge-Kutta 法求解，读者可将两个结果在一个图中显示，并和精确值比较，看看哪个更“精确”？（相应的 Matlab 程序参见附录 2）．四、自己动手1. 求微分方程的通解．0sin 2')1(2=-+-x xy y x 2. 求微分方程的通解．x e y y y x sin 5'2''=+-3. 求微分方程组⎪⎪⎩⎪⎪⎨⎧=-+=++00y x dtdy y x dt dx在初始条件下的特解，并画出解函数的图形．0|,1|00====t t y x ()y f x =4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解)，求解区间为．利用画图来比较两种求解器之间的差异．[0,2]t ∈5. 用 Euler 折线法求解微分方程初值问题⎪⎩⎪⎨⎧=-=1)0(,12'32y y x y y 的数值解(步长h 取0.1)，求解范围为区间[0,2]．6. 用四阶 Runge-Kutta 法求解微分方程初值问题⎩⎨⎧=-=1)0(,cos 'y x e y y x 准、流思，合格员”学、的数值解(步长h 取0.1)，求解范围为区间[0,3]．四阶 Runge-Kutta 法的迭代公式为(Euler 折线法实为一阶 Runge-Kutta 法)：1,,2,1,0),()2,2()2,2(),()22(6,),(342312143211100-=⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎨⎧++=++=++==++++=+==++n k hL y h x f L L h y h x f L L h y h x f L y x f L L L L L hy y h x x x y y k k kk k k k k k k k k 相应的 Matlab 程序参见附录 2．试用该方法求解第5题中的初值问题．7. 用 ode45 方法求上述第 6 题的常微分方程初值问题的数值解(近似解)，从而利用画图来比较两者间的差异．五、附录附录 1：(fulu1.m)clearf=sym('y+2*x/y^2');a=0;b=2;h=0.4;n=(b-a)/h+1;x=0;y=1;szj=[x,y];for i=1:n-1y=y+h*subs(f,{'x','y'},{x,y});x=x+h;szj=[szj;x,y];end szjplot(szj(:,1),szj(:,2))附录 2：(fulu2.m)clearf=sym('y-exp(x)*cos(x)');a=0;b=3;h=0.1;n=(b-a)/h+1;x=0;y=1;szj=[x,y];for i=1:n-1l1=subs(f,{'x','y'},{x,y});l2=subs(f,{'x','y'},{x+h/2,y+l1*h/2});l3=subs(f,{'x','y'},{x+h/2,y+l2*h/2});l4=subs(f,{'x','y'},{x+h,y+l3*h});y=y+h*(l1+2*l2+2*l3+l4)/6;x=x+h;szj=[szj;x,y];endszjplot(szj(:,1),szj(:,2))。
matlab解四阶偏微分
matlab解四阶偏微分在Matlab中，可以使用偏微分方程来解决四阶偏微分方程。
在本文中，我们将介绍四阶偏微分方程的一般形式、数值解法和一些相关的参考材料。
四阶偏微分方程的一般形式为：D^4u(x,y) + a*D^2u(x,y) + bu(x,y) = f(x,y)其中，D^4表示四阶空间导数算子，a和b是常数项，u(x,y)是要求解的未知函数，f(x,y)是已知的函数。
在Matlab中，可以使用偏微分方程工具箱（Partial Differential Equation Toolbox）来求解这个方程。
偏微分方程工具箱提供了多种数值方法来解决偏微分方程，包括有限差分法、有限元法、伽辽金法等。
有限差分法是最常用的数值方法之一，它将偏微分方程转化为一组有限差分方程，然后使用迭代方法求解这组方程。
有限差分法的基本思想是将求解区域离散化为网格，然后在网格节点上近似表示未知函数和导数。
通过在节点上构造差分方程，可以得到一个线性方程组，然后使用迭代方法求解这个方程组。
除了有限差分法，偏微分方程工具箱还提供了其他数值方法。
例如，有限元法将求解区域划分为多个小区域，然后在每个小区域内近似表示未知函数。
通过构造一组局部方程和边界条件，可以得到一个大型的线性方程组，然后使用迭代方法求解。
伽辽金法是一种通过变分原理求解偏微分方程的方法，它通过选取一个合适的试验函数，将偏微分方程转化为一组变分方程，然后通过极小化泛函来求解。
在Matlab中，偏微分方程工具箱提供了丰富的函数和工具来求解四阶偏微分方程。
例如，可以使用pdepe函数来求解带有边界条件的四阶偏微分方程，可以使用pdenonlin函数来求解非线性四阶偏微分方程。
此外，偏微分方程工具箱还提供了可视化工具和后处理函数，可以将求解结果可视化并进行进一步的分析。
除了Matlab自带的偏微分方程工具箱，还有一些其他的参考材料可以帮助理解和求解四阶偏微分方程。
例如，《Partial Differential Equations for Scientists and Engineers》是一本经典的偏微分方程教材，介绍了偏微分方程的基本理论和求解方法。
matlab解偏微分方程
matlab解偏微分方程Matlab是一种非常强大的数学计算工具，它可以用于解决各种数学问题。
在本文中，我们将学习如何使用Matlab解偏微分方程。
偏微分方程是一类包含未知函数的偏导数的方程。
通常，解偏微分方程是困难的，需要使用复杂的数学方法。
然而，Matlab可以大大简化这个过程。
在Matlab中，我们可以使用pdepe函数来解偏微分方程。
pdepe函数采用一个偏微分方程的系统，并返回一个包含解的向量的矩阵。
下面是一个解二维扩散方程的示例程序：%定义二维扩散方程 function [c,f,s] = diffusionpde(x,t,u,DuDx)c = 1; %系数f = DuDx; %带有时间和空间导数的项s = 0; %不带导数的项end%定义边界条件（例）function [pl,ql,pr,qr] =diffusionbc(xl,ul,xr,ur,t)pl = 0; ql = 1; %左边界（u=0）pr = 0; qr = 1; %右边界（u=0）end%定义初始条件（例）function u0 = diffusionic(x)u0 = sin(pi*x); %sin(pi*x)是初始条件方程end%主程序x = linspace(0,1,50); %空间网格t = linspace(0,1,10); %时间网格sol =pdepe(0,@diffusionpde,@diffusionic,@diffusionbc,x,t );u = sol(:,:,1); %提取第一个解%绘制解surfc(x,t,u)xlabel('位置')ylabel('时间')title('二维扩散方程的解')从上述程序中，我们可以看到pdepe的使用方法。
在主程序中，我们选择了空间和时间网格，然后定义了偏微分方程、初始条件和边界条件的函数。
最后，我们调用pdepe函数，并将解存储在变量sol中。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
1 / 16 [原创]偏微分方程数值解法的MATLAB源码【更新完毕】 说明：由于偏微分的程序都比较长，比其他的算法稍复杂一些，所以另开一贴，专门上传偏微分的程序 谢谢大家的支持！
其他的数值算法见： ..//Announce/Announce.asp?BoardID=209&id=8245004
1、古典显式格式求解抛物型偏微分方程（一维热传导方程）
function [U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C) %古典显式格式求解抛物型偏微分方程 %[U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C) % %方程：u_t=C*u_xx 0 <= x <= uX,0 <= t <= uT %初值条件：u(x,0)=phi(x) %边值条件：u(0,t)=psi1(t), u(uX,t)=psi2(t) % %输出参数：U -解矩阵，第一行表示初值，第一列和最后一列表示边值，第二行表示第2层…… % x -空间变量 % t -时间变量 %输入参数：uX -空间变量x的取值上限 % uT -时间变量t的取值上限 % phi -初值条件，定义为内联函数 % psi1 -边值条件，定义为内联函数 % psi2 -边值条件，定义为内联函数 % M -沿x轴的等分区间数 % N -沿t轴的等分区间数 % C -系数，默认情况下C=1 % %应用举例： %uX=1;uT=0.2;M=15;N=100;C=1; %phi=inline('sin(pi*x)');psi1=inline('0');psi2=inline('0'); %[U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C);
%设置参数C的默认值 if nargin==7 C=1; end
%计算步长 dx=uX/M;%x的步长 dt=uT/N;%t的步长 2 / 16
x=(0:M)*dx; t=(0:N)*dt;
r=C*dt/dx/dx;%步长比 r1=1-2*r;
if r > 0.5 disp('r > 0.5,不稳定') end
%计算初值和边值 U=zeros(M+1,N+1); for i=1:M+1 U(i,1)=phi(x(i)); end for j=1:N+1 U(1,j)=psi1(t(j)); U(M+1,j)=psi2(t(j)); end
%逐层求解 for j=1:N for i=2:M U(i,j+1)=r*U(i-1,j)+r1*U(i,j)+r*U(i+1,j); end end
U=U'; %作出图形 mesh(x,t,U); title('古典显式格式，一维热传导方程的解的图像') xlabel('空间变量 x') ylabel('时间变量 t') zlabel('一维热传导方程的解 U')
return;
古典显式格式不稳定情况 3 / 16
古典显式格式稳定情况 2、古典隐式格式求解抛物型偏微分方程（一维热传导方程） function [U x t]=PDEParabolicClassicalImplicit(uX,uT,phi,psi1,psi2,M,N,C) %古典隐式格式求解抛物型偏微分方程 %[U x t]=PDEParabolicClassicalImplicit(uX,uT,phi,psi1,psi2,M,N,C) % %方程：u_t=C*u_xx 0 <= x <= uX,0 <= t <= uT %初值条件：u(x,0)=phi(x) %边值条件：u(0,t)=psi1(t), u(uX,t)=psi2(t) % %输出参数：U -解矩阵，第一行表示初值，第一列和最后一列表示边值，第二行表示第2层…… % x -空间变量 % t -时间变量 4 / 16
%输入参数：uX -空间变量x的取值上限 % uT -时间变量t的取值上限 % phi -初值条件，定义为内联函数 % psi1 -边值条件，定义为内联函数 % psi2 -边值条件，定义为内联函数 % M -沿x轴的等分区间数 % N -沿t轴的等分区间数 % C -系数，默认情况下C=1 % %应用举例： %uX=1;uT=0.2;M=50;N=50;C=1; %phi=inline('sin(pi*x)');psi1=inline('0');psi2=inline('0'); %[U x t]=PDEParabolicClassicalImplicit(uX,uT,phi,psi1,psi2,M,N,C);
%设置参数C的默认值 if nargin==7 C=1; end
%计算步长 dx=uX/M;%x的步长 dt=uT/N;%t的步长
x=(0:M)*dx; t=(0:N)*dt;
r=C*dt/dx/dx;%步长比 Diag=zeros(1,M-1);%矩阵的对角线元素 Low=zeros(1,M-2);%矩阵的下对角线元素 Up=zeros(1,M-2);%矩阵的上对角线元素 for i=1:M-2 Diag(i)=1+2*r; Low(i)=-r; Up(i)=-r; end Diag(M-1)=1+2*r;
%计算初值和边值 U=zeros(M+1,N+1); for i=1:M+1 U(i,1)=phi(x(i)); end 5 / 16
for j=1:N+1 U(1,j)=psi1(t(j)); U(M+1,j)=psi2(t(j)); end
%逐层求解，需要使用追赶法（调用函数EqtsForwardAndBackward） for j=1:N b1=zeros(M-1,1); b1(1)=r*U(1,j+1); b1(M-1)=r*U(M+1,j+1); b=U(2:M,j)+b1; U(2:M,j+1)=EqtsForwardAndBackward(Low,Diag,Up,b); end
U=U'; %作出图形 mesh(x,t,U); title('古典隐式格式，一维热传导方程的解的图像') xlabel('空间变量 x') ylabel('时间变量 t') zlabel('一维热传导方程的解 U')
return;
此算法需要使用追赶法求解三对角线性方程组，这个算法在上一篇帖子中已经给出，为了方便，再给出来 追赶法解三对角线性方程组
function x=EqtsForwardAndBackward(L,D,U,b) %追赶法求解三对角线性方程组Ax=b %x=EqtsForwardAndBackward(L,D,U,b) %x:三对角线性方程组的解 %L:三对角矩阵的下对角线，行向量 %D:三对角矩阵的对角线，行向量 %U:三对角矩阵的上对角线，行向量 %b:线性方程组Ax=b中的b，列向量 % %应用举例: %L=[-1 -2 -3];D=[2 3 4 5];U=[-1 -2 -3];b=[6 1 -2 1]'; 6 / 16
%x=EqtsForwardAndBackward(L,D,U,b) %检查参数的输入是否正确 n=length(D);m=length(b); n1=length(L);n2=length(U); if n-n1 ~= 1 || n-n2 ~= 1 || n ~= m disp('输入参数有误！') x=' '; return; end
%追的过程 for i=2:n L(i-1)=L(i-1)/D(i-1); D(i)=D(i)-L(i-1)*U(i-1); end x=zeros(n,1); x(1)=b(1); for i=2:n x(i)=b(i)-L(i-1)*x(i-1); end
%赶的过程 x(n)=x(n)/D(n); for i=n-1:-1:1 x(i)=(x(i)-U(i)*x(i+1))/D(i); end
return;
古典隐式格式 7 / 16
在以后的程序中，我们都取C=1，不再作为一个输入参数处理 3、Crank-Nicolson隐式格式求解抛物型偏微分方程 需要调用追赶法的程序
function [U x t]=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N) %Crank-Nicolson隐式格式求解抛物型偏微分方程 %[U x t]=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N) % %方程：u_t=u_xx 0 <= x <= uX,0 <= t <= uT %初值条件：u(x,0)=phi(x) %边值条件：u(0,t)=psi1(t), u(uX,t)=psi2(t) % %输出参数：U -解矩阵，第一行表示初值，第一列和最后一列表示边值，第二行表示第2层…… % x -空间变量 % t -时间变量 %输入参数：uX -空间变量x的取值上限 % uT -时间变量t的取值上限 % phi -初值条件，定义为内联函数 % psi1 -边值条件，定义为内联函数 % psi2 -边值条件，定义为内联函数 % M -沿x轴的等分区间数 % N -沿t轴的等分区间数 % %应用举例： %uX=1;uT=0.2;M=50;N=50; %phi=inline('sin(pi*x)');psi1=inline('0');psi2=inline('0'); %[U x t]=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N);