数学建模 两点边值问题的两种数值解法及代码
重要：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];
数学建模介绍
数学建模介绍1.1 数学模型及其分类数学建模作为用数学方法解决问题的第一步，它与数学本身有着同样悠久的历史。
一个羊倌看着他的羊群进入羊圈，为了确信他的羊没有丢失，他在每只羊进入羊圈时，则在旁边放一颗小石子，如果每天羊全部入圈而他那堆小石子刚好全部放完，则表示他的羊和以前一样多。
究竟羊倌数的是石子还是羊，那是毫无区别的，因为羊的数目同石子的数目彼此相等。
这实际上就使石子与羊“联系”起来，建立了一个使石子与羊一一对应的数学模型。
（1）什么是数学模型人们在认识研究现实世界里的客观对象时，常常不是直接面对那个对象的原形，有些是不方便，有些甚至是不可能直接面对原形，因此，常常设计、构造它的各种各样的模型。
如各式各样的玩具模型、展览厅里的三峡大坝模型、化学上的分子结构模型等。
这些模型都是人们为了一定目的，对客观事物的某一部分进行简化、抽象、提炼出来的原形替代物，集中反映了原形中人们需要的那一部分特征，因而有利于人们对客观对象的认识。
数学模型也是反映客观对象特征的，只不过它刻画的是事物在数量方面的特征或数学结构及其变化规律。
数学模型是人们为了认识客观对象在数量方面的特征、定量地分析对象的内在规律、用数学的语言和符号去近似地刻画要研究的那一部分现象时，所得到的一个数学表述。
建立数学模型的过程称为数学建模。
(2) 数学模型的重要作用进入20世纪以来，数学以空前的广度和深度向一切领域渗透，作为数学的应用，数学建模也越来越受到人们的重视。
在一般工程技术领域，数学模型仍是工程技术人员定量研究有关工程技术问题的重要工具；而随着数学与其他学科领域诸如经济、人口、生态、地质等所谓非物理领域的渗透，一些交叉学科如计量经济学、人口控制论、数学生态学、数学地质学等应运而生；计算机的发展给数学及作为数学应用的数学建模带来了前所未有的机遇和挑战。
计算机改变了人类的生活方式、思考方式和研究方式，极大地提高了人们的计算能力、搜索和分析海量数据和信息的能力。
数学建模中的图论方法
数学建模中的图论方法一、引言我们知道，数学建模竞赛中有问题A和问题B。
一般而言，问题A是连续系统中的问题，问题B是离散系统中的问题。
由于我们在大学数学教育内容中，连续系统方面的知识的比例较大，而离散数学比例较小。
因此很多人有这样的感觉，A题入手快，而B题不好下手。
另外，在有限元素的离散系统中，相应的数学模型又可以划分为两类，一类是存在有效算法的所谓P类问题，即多项式时间内可以解决的问题。
但是这类问题在MCM中非常少见，事实上，由于竞赛是开卷的，参考相关文献，使用现成的算法解决一个P类问题，不能显示参赛者的建模及解决实际问题能力之大小；还有一类所谓的NP问题，这种问题每一个都尚未建立有效的算法，也许真的就不可能有有效算法来解决。
命题往往以这种NPC问题为数学背景，找一个具体的实际模型来考验参赛者。
这样增加了建立数学模型的难度。
但是这也并不是说无法求解。
一般来说，由于问题是具体的实例，我们可以找到特殊的解法，或者可以给出一个近似解。
图论作为离散数学的一个重要分支，在工程技术、自然科学和经济管理中的许多方面都能提供有力的数学模型来解决实际问题，所以吸引了很多研究人员去研究图论中的方法和算法。
应该说，我们对图论中的经典例子或多或少还是有一些了解的，比如，哥尼斯堡七桥问题、中国邮递员问题、四色定理等等。
图论方法已经成为数学模型中的重要方法。
许多难题由于归结为图论问题被巧妙地解决。
而且，从历年的数学建模竞赛看，出现图论模型的频率极大，比如：AMCM90B－扫雪问题；AMCM91B－寻找最优Steiner树；AMCM92B－紧急修复系统的研制(最小生成树)AMCM94B－计算机传输数据的最小时间(边染色问题)CMCM93B－足球队排名(特征向量法)CMCM94B－锁具装箱问题(最大独立顶点集、最小覆盖等用来证明最优性)CMCM98B－灾情巡视路线(最优回路)等等。
这里面都直接或是间接用到图论方面的知识。
要说明的是，这里图论只是解决问题的一种方法，而不是唯一的方法。
数学建模(7)
a1 j a Aj = 2j ... a mj
a11a12 ...a1n a a ...a A = ( A1 , A2 ,..., An ) = 21 22 2 n ... a a ...a m 1 m 2 mn
优化模型
优化模型
3、单变量优化的Matlab实现 、单变量优化的 实现 fminbnd 返回固定区间内单变量函数的最小值 x = fminbnd(fun,x1,x2) 返回区间{x1， 上 参数描述的标量函数的最小 返回区间 ，x2}上fun参数描述的标量函数的最小 值x x = fminbnd(fun,x1,x2,options) 用options参数指定的优化参数进行最小化 参数指定的优化参数进行最小化 x = fminbnd(fun,x1,x2,options,P1,P2,...) 提供另外参数P 传输给目标函数fun。 提供另外参数 1,P2等,传输给目标函数 。如果没有 传输给目标函数 设置options选项，则令 选项， 设置 选项 则令options=[]
mq (α ) = aα + bα + c
2
其中步长极值为： 其中步长极值为：
−b α = 2a
∗
优化模型
然后只要利用三个梯度或函数方程组就可以确定系 数a和b，从而可以确定 。得到该值以后，进行搜索区 和 ，从而可以确定α*。得到该值以后， 间的收缩。在缩短的新区间中， 间的收缩。在缩短的新区间中，重新安排三点求出下一 次的近似极小点α*，如此迭代下去， 次的近似极小点 ，如此迭代下去，直到满足终止准则 为止。其迭代公式为： 为止。其迭代公式为：
（4）线性规划问题的标准形式是： ）线性规划问题的标准形式是： n max z = c1 x1 + c2 x 2 + ⋅ ⋅ ⋅ + cn x n max z = ∑ c j x j a x + a x + ⋅⋅⋅ + a x = b j =1 n 1n n 1 11 1 12 2 或 ∑ aij x j = bi , i = 1,2,..., m a 21 x1 + a 22 x 2 + ⋅ ⋅ ⋅ + a 2 n x n = b2 j =1 x j ≥ 0, j = 1,2,..., n ⋅⋅⋅⋅⋅⋅ a m 1 x1 + a m 2 x 2 + ⋅ ⋅ ⋅ + a mn x n = bm x1 , x 2 ,⋅ ⋅ ⋅, x n ≥ 0 矩阵形式为： 矩阵形式为：
数学建模第八讲：偏微分方程数值解
2 (t )
其中：u
t
0
(
x
),
u t
t0
(x)
为初值条件
u x0 1 (t ), u xt 2 (t ) 为边值条件
当该波动方程只提供初值条件时，称此方程为波动方程的初值问题，二
。 者均提供时称为波动方程的混合问题
5.3.1 波动方程求解
t
t
x 0 a)初值问题
x
0
l
b)混合问题
对于初值问题，是已知t=0时，u与u 依赖于x的函数形式，求解不同位置， t
un1 i , j,k
t 2 t nt , xix , y jy,zkz
( t )2
2u x 2
t nt , xix , y jy,zkz
un i1, j,k
2uin, j,k (x)2
un i1, j ,k
2u y2
t nt , xix , y jy,zkz
un i , j1,k
2uin, j,k (y)2
21
A11 I
A
I
A22 I
I AN 2 ,N 2 I
I
R( N 1)2 ( N 1)2
AN 1,N 1
其中
4 1
f ( x, t)
u
t
0
(
x
),
u t
t0
( x)
u
x0
1(t), u xl
2(t)
uin
un1 i
τn
xi
x
un1 i
方程离散化
un1 i
2uin
un1 i
(t )2
a2
un i1
2uin
python数值求解微分方程
一、介绍1. Python是一种高级编程语言，具有简单易学、功能强大、开源免费等优点。
在科学计算、数据分析、机器学习等领域被广泛应用。
2. 微分方程是数学中的一个重要分支，描述了变化的规律。
微分方程求解是许多科学工程问题的基础和关键。
二、微分方程的数值求解1. 微分方程的数值求解是通过数值方法近似计算微分方程的解。
常用的数值方法包括欧拉方法、隐式欧拉方法、龙格-库塔方法等。
2. 对于一阶微分方程y'=f(x,y)，欧拉方法的迭代公式为y_n+1=y_n+hf(x_n,y_n)，其中h为步长。
3. 隐式欧拉方法则采用y_n+1=y_n+hf(x_n+1,y_n+1)，需要通过迭代求解非线性方程。
4. 龙格-库塔方法是一种更精确的数值方法，通常采用四阶龙格-库塔方法。
三、使用Python进行微分方程数值求解1. Python具有丰富的数值计算库，如NumPy、SciPy等，能方便地进行微分方程的数值求解。
2. 使用SciPy库中的odeint函数可以方便地对常微分方程进行数值求解。
3. 以一阶常微分方程y'=2x+y，y(0)=1为例，可以使用Python实现如下：import numpy as npfrom scipy.integrate import odeintdef model(y, x):return 2*x + yx = np.linspace(0, 5, 100)y = odeint(model, 1, x)4. 这段代码首先定义了微分方程的模型函数model，并使用odeint 函数进行数值求解。
最终得到的y即为微分方程的近似解。
四、实例分析1. 以一个具体的物理问题为例，考虑一个简谐振动系统，其运动方程为mx''+kx=0。
假设m=1，k=1，初始条件为x(0)=0，x'(0)=1。
2. 可以将这个二阶微分方程转化为一对一阶微分方程，令y1=x，y2=x'，则转化得到y1'=y2，y2'=-y1。
915225-MATLAB程序设计与应用-第8章 MATLAB方程数值求解__源程序
第8章MATLAB方程数值求解例8-1用直接解法求解下列线性方程组。
程序如下：A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';x=A\b例8-2用LU分解求解例8-1中的线性方程组。
程序如下：A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';[L,U]=lu(A);x=U\(L\b)例8-3 用QR分解求解例8-1中的线性方程组。
程序如下：A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';[Q,R]=qr(A);x=R\(Q\b)例8-4 用Cholesky分解求解例8-1中的线性方程组。
命令如下：>> A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; >> b=[13,-9,6,0]';>> R=chol(A)Jacobi迭代法的MA TLAB函数文件jacobi.m如下：function [y,n]=jacobi(A,b,x0,ep)if nargin==3ep=1.0e-6;elseif nargin<3errorreturnendD=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵B=D\(L+U);f=D\b;y=B*x0+f;n=1; %迭代次数while norm(y-x0)>=epx0=y;y=B*x0+f;n=n+1;end例8-5 用Jacobi迭代法求解下列线性方程组。
设迭代初值为0，迭代精度为10-6。
在程序中调用函数文件jacobi.m，程序如下：A=[10,-1,0;-1,10,-2;0,-2,10];b=[9,7,6]';[x,n]=jacobi(A,b,[0,0,0]',1.0e-6)Gauss-Serdel迭代法的MA TLAB函数文件gauseidel.m如下：function [y,n]=gauseidel(A,b,x0,ep)if nargin==3ep=1.0e-6;elseif nargin<3errorreturnendD=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵G=(D-L)\U;f=(D-L)\b;y=G*x0+f;n=1; %迭代次数while norm(y-x0)>=epx0=y;y=G*x0+f;n=n+1;end例8-6 用Gauss-Serdel迭代法求解例8-5中的线性方程组。
半解析法 与半解析数值法的缩写
半解析法与半解析数值法的缩写半解析法和半解析数值法是数值计算中常用的两种方法，它们在求解一些特定类型的数学问题时具有一定的优势。
下面将对半解析法和半解析数值法的定义、原理和应用进行详细说明。
半解析法（Semi-Analytical Method）是指通过对问题进行分析和推导，得到问题的一种解析形式，但仍然需要使用数值计算方法来求解的方法。
与传统的纯解析法相比，半解析法可以更好地处理一些复杂的问题，因为它可以通过分析和推导将问题转化为更简单的形式。
半解析法的应用十分广泛，特别是在物理学、工程学和数学建模等领域。
例如，在电磁场计算中，经常需要通过求解电磁场的边值问题来得到电磁场的分布情况。
传统的纯解析法很难处理复杂边界条件下的电磁场问题，而半解析法可以通过数学推导得到电磁场的一种解析形式，并进一步使用数值计算方法得到精确的数值解。
半解析法的原理基于分析和推导，通常可以采用比较精确的数学工具，如微分方程、积分变换和级数展开等。
通过对问题进行适当的简化和近似，可以得到问题的解析形式。
然后，可以使用数值计算方法对解析形式进行进一步处理，求取数值解。
与半解析法相对应的是半解析数值法（Semi-Analytical Numerical Method）。
半解析数值法是指将数值计算方法与解析方法相结合的一种数值计算方法。
它通过分析和推导的方法，得到问题的解析形式，并且使用数值计算方法对解析形式进行进一步处理，求取数值解。
半解析数值法的优势在于既保留了解析法的精确性，又克服了纯解析法无法处理复杂问题的困难。
半解析数值法在求解偏微分方程、积分方程和边值问题等数学问题时具有很大的优势。
半解析数值法的应用范围广泛，包括流体力学、结构力学、电磁场计算、热传导等领域。
例如，在流体力学中，常常需要求解Navier-Stokes方程来描述流体的运动行为。
传统的纯解析法很难处理复杂的流动问题，而半解析数值法可以通过求解Navier-Stokes方程的一种解析形式，并通过数值计算方法对解析形式进行进一步处理，得到精确的数值解。
高等应用数学问题的MATLAB求解课件
随着科技的发展，数学在各个领域的应用越来越广泛，而高等应用数学作为数学的一个重要分支，其 重要性不言而喻。为了更好地理解和应用高等应用数学，我们需要掌握如何使用工具软件进行求解。
目标
本课程旨在帮助学生掌握使用Matlab软件解决高等应用数学问题的方法，提高分析和解决实际问题的 能力。
Matlab软件简介
• 详细描述：在Matlab中，可以通过设置不同的参数来控制算法的收敛速度、 精度等。例如，fmincon函数中的'Algorithm'参数可以设置为'sqp'、'trustregion'等，以适应不同的问题规模和复杂度。
矩阵特征值问题
• 总结词：矩阵特征值问题是一类重要的数学问题，它涉及到矩阵的特征值和特 征向量。
• 详细描述：在Matlab中，可以使用eig、eigs等函数来求解矩阵的特征值和特 征向量。这些函数可以处理各种类型的矩阵，包括实数矩阵、复数矩阵等。通 过计算矩阵的特征值和特征向量，可以解决许多实际问题，如振动分析、控制 系统设计等。
• 总结词：求解矩阵特征值问题时，需要注意数值稳定性问题，以避免计算误差 和数值不稳定性。
05
Matlab在数学建模中的应用
建模基础
变量与数据
确定问题中的变量和数据，为建模提供基础。
数学关系建立
根据问题背景和实际需求，建立数学关系式， 描述变量间的关系。
模型简化与求解
对建立的数学模型进行简化，并使用 Matlab进行求解。
建模实例分析
01
线性回归模型
02
非线性拟合模型
03
微分方程求解
复杂系统模拟
离散事件模拟
01
利用Matlab进行离散事件模拟，适用于模拟具有离散时间或状
数学建模图论讲
第2页1 /共86页
2024年8月3日
数学建模-图论
一、图的基本概念
如果图的二顶点间有边相连，则称此顶点相邻，每一对顶点
都相邻的图称为完全图，否则称为非完全图，完全图记为 K V 。
若V (G) X Y, X Y , X Y 0 ，且 X 中 无相邻的顶点对，Y 中亦然，则称图 G 为二分图.
第1行 1 A1i 第i行 1
11，A1i 2
2 2
22，A1i3
4 4
4 4
其中i=2,3,4,5,显然y1=1+(4+4+4+4-1) 4=61. 同理，计算y2时应考虑槽高只有2,21,23,24,25,
26时的情形，类似计算可得 y2=1+(4+4+4+4-1)×5=76．
于是，s=61×2+76×4=426，x=6306426=5880.
计算y1可分别考虑槽高只有1，12，13，14，15的 情形．若只有1，这样的锁具效只有1个， 若只有1和i(i=2,3,4,5)，这样的锁具数=G中以1和i为 顶点，长度为3的道路数，此数可通过A的子矩阵A1i计 算得到．
第18页/共86页
数学建模-图论
二、图的矩阵表示（应用实例解法分析）
事实上，因为
间最短的路线。定义T*T=(t(2)ij),
3
4
t(2)ij=min{min1<=k<=5{tik+tkj},tij}, t(2)ij表示 从站点i到站点j的至多换乘一次的最短时间。
5
第22页/共86页
数学建模-图论
二、图的矩阵表示（应用实例及解法分析）
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
常微分方程组两点边值问题的数值解法
3)1(1)0(04y
y
yy
可化为微分方程组3)1(1)0(41221yyyyyy
方法一：配置法
Matlab程序：
function bvcollation
clc
solinit = bvpinit(linspace(0,1,20),[100 600]);%
sol = bvp4c(@twoode,@twobc,solinit);
x = linspace(0,1,20);
y = deval(sol,x);
y'
plot(x,y(1,:),x,y(2,:));
end
%微分方程组
function dydx = twoode(x,y)
dydx = [ y(2)
4*y(1)];
end
%边值条件
function res = twobc(ya,yb)
res = [ ya(1)-1
yb(1)-3];
end
运行结果：
1.0000 -0.4203
0.9834 -0.2117
0.9777 -0.0055
0.9828 0.2007
0.9988 0.4091
1.0259 0.6220
1.0644 0.8419
1.1147 1.0710
1.1774 1.3121
1.2531 1.5677
1.3427 1.8407
1.4472 2.1341
1.5678 2.4512
1.7057 2.7954
1.8626 3.1707
2.0401 3.5811
2.2402 4.0313
2.4652 4.5261
2.7175 5.0712
3.0000 5.6724
方法二：打靶法
程序：积分用ode45,搜索用二分法
function shoot_first_try
clc
clear all
a=0;
b=1;
a2=-1;
b2=0;
s=0;
tol=1.0*10e-6;
h1=Fun(a2)
h2=Fun(b2)
s=bisect(@Fun,a2,b2,tol)
[t1,y1]=ode45(@f,[a,b],[1,s])
plot(t1,y1(:,1),t1,y1(:,2))
function xc=bisect(f,a,b,tol)
if sign(f(a))*sign(f(b))>=0
error('f(a)*f(b)<0 not satisfied!')
end
fa=f(a);
fb=f(b);
k=0;
while (b-a)/2>tol
c=(a+b)/2;
fc=f(c);
if fc==0
break
end
if sign(f(a))*sign(f(b))<0
b=c;fb=fc;
else
a=c;fa=fc;
end
end
xc=(a+b)/2;
function z=Fun(s)
a=0;b=1;yb=3;
h=linspace(a,b,20);
y0=[1;s];
[t,y]=ode45(@f,h,y0);
z=y(end,1)-yb;
function ydot=f(t,y)
%ydot=[0;0];
ydot(1)=y(2);
ydot(2)=4*y(1);
ydot=[ydot(1);ydot(2)];
运行结果：
1.0000 -0.5000
0.9969 -0.4749
0.9940 -0.4499
0.9913 -0.4250
0.9887 -0.4001
0.9799 -0.3017
0.9736 -0.2041
0.9697 -0.1069
0.9683 -0.0100
0.9692 0.0868
0.9726 0.1839
0.9784 0.2814
0.9867 0.3797
0.9974 0.4788
1.0106 0.5792
1.0264 0.6811
1.0447 0.7846
1.0656 0.8901
1.0892 0.9978
1.1155 1.1080
1.1446 1.2210
1.1766 1.3370
1.2115 1.4564
1.2495 1.5794
1.2905 1.7064
1.3348 1.8377
1.3825 1.9735
1.4335 2.1143
1.4882 2.2603
1.5466 2.4120
1.6089 2.5698
1.6751 2.7339
1.7456 2.9049
1.8204 3.0832
1.8998 3.2692
1.9840 3.4633
2.0731 3.6661
2.1674 3.8781
2.2671 4.0998
2.3724 4.3317
2.4837 4.5745
2.5711 4.7637
2.6621 4.9596
2.7569 5.1625
2.8555 5.3726