bvp4c--语法

 高质量学习资源免费获取,专注但不限于【Linux】【C/C++/Qt】【FPGA】【数据结构与算法】, 根据多年技术经验纯【原创】,纯【干货】,分享【技术心得】,力求【授人以鱼,更授人以渔】。

bvp4c--语法  

1. bvp4c

sol = bvp4c(odefun,bcfun,solinit)

sol = bvp4c(odefun,bcfun,solinit,options)

sol = bvp4c(odefun,bcfun,solinit,options,p1,p2...)

sol返回如下值

sol.x  bvp4c选择的网格

sol.y   在网格点sol.x的 y(x)的接近值

sol.yp 在网格点sol.x的 y'(x)的接近值

sol.parameters 未知参数的值

sol.solver   'bvp4c'

2.odefun 为微分方程

dydx = odefun(x,y)
dydx = odefun(x,y,p1,p2,...)
dydx = odefun(x,y,parameters)
dydx = odefun(x,y,parameters,p1,p2,...)

parameters 为未知参数向量,p1和p2...为已知参数向量

dydx为列向量

3. bcfun 计算边界条件的残留值bc(  y(a),y(b)  ),输出列向量

res = bcfun(ya,yb)
res = bcfun(ya,yb,p1,p2,...)
res = bcfun(ya,yb,parameters)
res = bcfun(ya,yb,parameters,p1,p2,...)

4. solinit

solinit的结构如下:

x  初始网格的节点(相当于自变量的节点),边界条件满足a=solinit.x(1) 和 b=solint.x(end)

y  解的初始假设,solinit.y(:,i) 为自变量solinit.x(i)对应的解的估计值

parameters 可选,未知参数的初始估计向量,对于存在未知参数的情况,必须提供

sol.y(1,1).............. sol.y(1,i)............. sol.y(1,end)

...

sol.y(m,1).............. sol.y(m,i)............. sol.y(m,end)

sol.x(1)   ................sol.x(i)  .................sol.x(end)

solinit的赋值语句为bvpinit,语法如下

solinit = bvpinit(x,v)
solinit = bvpinit(x,v,parameters)
solinit = bvpinit(sol,[anew bnew])
solinit = bvpinit(sol,[anew bnew],parameters)

x为初始网格向量,对于边界区间[a,b]需满足 x(1)=a, x(end)=b, 一般情况下可用 x=linspace(a,b,10)

v为解的估计,可以是一个向量,也可以是一个函数

v为向量的情况,v(i)代表的就是 y(i,:),也就是不管自变量为多少y(i)都是这个值,向量的维数等于因变量个数

v为函数的情况,对于给定一个网格,函数必须返回一个向量y(1),y(2),y(3)...y(m)

y=guess(x) 代表的是 y(:,j)=guess( x(j) )

solinit = bvpinit(sol,[anew bnew])

由[a,b]上的解sol,得到[anew,bnew]的解的初始估计,so either anew <= a < b <= bnew or anew >= a > b >= bnew

 5.deval 评估 自变量=xint时候的解

sxinit= deval(sol.xint)

6. options的设置 bvpset语句

options = bvpset('name1',value1,'name2',value2,...)
options = bvpset(oldopts'name1',value1,...)
options = bvpset(oldopts,newopts)

name               value

RelTol             正标量  默认为1e-3  相对精度

AbsTol            正标量 默认为1e-6  绝对精度

Vectorized      on或者off  将ode函数F([x1 x2 ...],[y1 y2 ...]),写成[F(x1,y1) F(x2,y2) ...]. 的形式,能减少函数评估的次数,降低运行时间

SingularTerm   矩阵    方程的奇异项, 如 y'=Sy/x + f(x,y,p)  x=[0,b] ,则设置常数矩阵S

FJacobian       函数、矩阵或元胞数组  提供odefun的解析偏导数 对于y'=f(x,y),则提供f对y的偏导数,若偏导数都为常数项,则用元胞数组的形式给出;若有未知参数,还需计算对未知参数的偏导数

BCJacobian    函数、元胞数组     提供边界函数bcfun的解析偏导数 bc(ya,yb)=0 则计算bc对ya,yb 的偏导数,若有未知参数p还需计算未知参数的偏导数:[DBCDYA,DBCDYB,DBCDP] = BCJAC(YA,YB,P)

Nmax              正整数    网格的最大值

Stats              on 或者off  显示计算过程的统计数据

 算法原理:

bvp4c is a finite difference code that implements the three-stage Lobatto IIIa formula. This is a collocation formula and the collocation polynomial provides a C1-continuous solution that is fourth order accurate uniformly in [a,b]. Mesh selection and error control are based on the residual of the continuous solution.

 1 % bvp4c_fun.m
 2 % 求解微分方程边值问题
 3 S=[ 0 0
 4     0 -2];
 5 options=bvpset('SingularTerm',S);
 6 % 得到微分方程的初场
 7 guess=[sqrt(3)/2;0];
 8 solinit=bvpinit(linspace(0,1,5),guess);
 9 sol=bvp4c(@emdenode,@emdenbc,solinit,options);
10 % 根据分析解计算
11 x=linspace(0,1);
12 truy=1./sqrt(1+(x.^2)/3);
13 % 绘制图形结果
14 plot(x,truy,'b-');
15 hold on;
16 plot(sol.x,sol.y(1,:),'ro');
17 legend('analytical','computed');
18 xlabel('x'),ylabel('y');
19 hold off;
View Code
原文地址:https://www.cnblogs.com/lihuidashen/p/3445001.html