【Matlab】三次B样条基函数插值求解泛函极值问题
---------------------------------------------------------------------------------------------
前言：博主是北理工的大三自动化人，分享学习经历，欢迎指导！
欢迎关注微信公众号：不断努力的霜木君
---------------------------------------------------------------------------------------------
前一阵子考完了科学与工程计算，但是作业还没写完Orz
看了熊老师布置的第四题，起初真的是一脸懵逼。
不是摸不着头脑，而是感觉这玩意也太复杂了吧！
天无绝人之路，在借鉴了骁哥的作业后，我有了做下去的动力。
事实上，我俩的方法一点都不一样，这里我记录一下我掉了好几根头发，按照最普通的方式搓出来的方法。
其中的关键点可能是……好像没啥关键点，就是生搓硬搓啊！
先说下题目，
很经典的一个求泛函极值的问题，求出Euler公式后，使用Galerkin方法最终化为
进而，考虑三次B样条函数
也就是求解
另外，本题的真解如下，用于进行误差比较
这样，问题基本上就解决了。现在只要构建方程组的矩阵形式来求解就可以了。
还有一个小细节要提一下，我这里为了计算方便，将20等分的区间（x标号为1~21）仅从i=3到i=19进行了B样条函数的构建。因为matlab里面序列是从1开始的。只不过影响不大，特别是区间更细之后就更没啥影响了。
当然，上述想法仅是我稍微偷了偷懒，如果真有问题的话（应该是有问题）可以公众号里发消息和我讨论。
最后得到的效果如下，蓝色的是真解，红色是近似解
1、区间等分20份
2、区间等分200份
3、区间等分1000份
可以看到，逼近效果随区间变细越来越好。
------------------------------------------------------------------------
接下来，就进入了编程过程。
具体不解释了，注释写的挺清楚。
没啥含金量，主要是积分采用了梯形法来近似。
可自行复制到matlab尝试，因为公众号上显示格式有问题
拜拜了大家，继续复习统计学了。
clear all
clc
N=1000;
h=1/N;
xlist=linspace(0,1,N+1);
ulist=zeros(N+1,1);
A=zeros(N-3,N-3);
LL=zeros(N-3,1);%构建矩阵求解方程组
for i=3:N-1
for j=3:N-1
A(i-2,j-2)=a(i,j,xlist);
end
LL(i-2)=L(i,xlist);
end
c=inv(A)*LL;
for i=3:N-1
for j=3:N-1
ulist(i-2)=ulist(i-2)+c(j-2)*phi(xlist,j,xlist(i-2),h);
end
end
res_real=sin(xlist)/sin(1)-xlist;
plot(xlist,ulist,'r');
hold on
plot(xlist,res_real,'b');
title(['区间',num2str(N),'等分的三次B样条插值结果']);
%三次B样条函数
function res = phi(xlist,i,x,h)
if x <= xlist(i-2)
res = 0;
elseif x>xlist(i-2) && x<=xlist(i-1)
res = (x-xlist(i-2))^3/(6*h^3);
elseif x>xlist(i-1) && x<=xlist(i)
res=(1/6)+(x-xlist(i-1))/(2*h)+(x-xlist(i-1))^2/(2*h^2)-(x-xlist(i-1))^3/(2*h^3);
elseif x>xlist(i) && x<=xlist(i+1)
res=(1/6)-(x-xlist(i+1))/(2*h)+(x-xlist(i+1))^2/(2*h^2)+(x-xlist(i+1))^3/(2*h^3);
elseif x>xlist(i+1) && x<=xlist(i+2)
res = -(x-xlist(i+2))^3/(6*h^3);
else
res = 0;
end
end
%三次B样条函数的微分
function res = dphi(xlist,i,x,h)
if x <= xlist(i-2)
res = 0;
elseif x>xlist(i-2) && x<=xlist(i-1)
res = (x-xlist(i-2))^2/(2*h^3);
elseif x>xlist(i-1) && x<=xlist(i)
res=(1)/(2*h)+(x-xlist(i-1))/(h^2)-3*(x-xlist(i-1))^2/(2*h^3);
elseif x>xlist(i) && x<=xlist(i+1)
res=-(1)/(2*h)+(x-xlist(i+1))/(h^2)+3*(x-xlist(i+1))^2/(2*h^3);
elseif x>xlist(i+1) && x<=xlist(i+2)
res = -(x-xlist(i+2))^2/(2*h^3);
else
res = 0;
end
end
%梯形积分公式对a(i,j) 0阶导部分积分
function I = Trapezoid(a,b,i, j,xlist)
dis=(b-a);
h=xlist(2)-xlist(1);
I=((dis)/2)*(phi(xlist,i,a,h)*phi(xlist,j,a,h)+phi(xlist,i,b,h)*phi(xlist,j,b,h));
end
%梯形积分公式对a(i,j) 1阶导部分积分
function I = dTrapezoid(a,b,i, j,xlist)
h=xlist(2)-xlist(1);
dis=(b-a);
I=((dis)/2)*(dphi(xlist,i,a,h)*dphi(xlist,j,a,h)+dphi(xlist,i,b,h)*dphi(xlist,j,b,h));
end
%梯形积分公式对L(i)积分
function I = LTrapezoid(a,b,i,xlist)
dis=(b-a);
h=xlist(2)-xlist(1);
I=((dis)/2)*(phi(xlist,i,a,h)*a+phi(xlist,i,b,h)*b);
end
%求解a(i,j)
function res = a(i,j,xlist)
if abs(i-j)>=4
res=0;
else
if i>j
temp=i;
i=j;
j=temp;
end
if (j-i) == 0
res=2*(Trapezoid(xlist(i-2),xlist(i-1),i,j,xlist)+Trapezoid(xlist(i-1),xlist(i),i, j,xlist))-2*(dTrapezoid(xlist(i-2),xlist(i-1),i, j,xlist)+dTrapezoid(xlist(i-1),xlist(i),i, j,xlist));
elseif (j-i)==1
res=Trapezoid(xlist(i-1),xlist(i),i,j,xlist)+Trapezoid(xlist(i),xlist(i+1),i, j,xlist)+Trapezoid(xlist(j),xlist(j+1),i,j,xlist)-(dTrapezoid(xlist(i-1),xlist(i),i, j,xlist)+dTrapezoid(xlist(i),xlist(i+1),i,j,xlist)+dTrapezoid(xlist(j),xlist(j+1),i, j,xlist));
elseif(j-i)==2
res=Trapezoid(xlist(i),xlist(i+1),i,j,xlist)*2-dTrapezoid(xlist(i),xlist(i+1),i, j,xlist)*2;
elseif(j-i)==3
res=Trapezoid(xlist(i+1),xlist(i+2),i,j,xlist)-dTrapezoid(xlist(i+1),xlist(i+2),i, j,xlist);
else
res=0;
end
end
end
%求解L(i)
function res=L(i,xlist)
res=LTrapezoid(xlist(i-2),xlist(i-1),i,xlist)+LTrapezoid(xlist(i-1),xlist(i),i,xlist)+LTrapezoid(xlist(i),xlist(i+1),i,xlist)+LTrapezoid(xlist(i+1),xlist(i+2),i,xlist);
res = res *(-1);
end
【Matlab】三次B样条基函数插值求解泛函极值问题相关推荐
- matlab怎么求三次微分,matlab课设三阶微分方程多种方法求解.doc
matlab课设三阶微分方程多种方法求解 目录 一.课程设计题目及意义 -------- 1 页 二.课程设计任务及要求 --------2 页 三.课程设计详细过程及结果 --------3至10页 ...
- matlab径向基函数插值,径向基函数（Radial Basis Function）插值
将RBF用于插值 标签(空格分隔):径向基函数插值 算法 RBF 曲面重构 当高维数据稀疏,需要预测一些数据,需要使用曲面重构的方法. 曲面重构一般可以分为: 插值 重构 曲面插值里我们一般使用径向基 ...
- 三次B样条曲线拟合算法
1 三次B样条曲线方程 B样条曲线分为近似拟合和插值拟合,所谓近似拟合就是不过特征点,而插值拟合就是通过特征点,但是插值拟合需要经过反算得到控制点再拟合出过特征点的B样条曲线方程.这里会一次介绍两种拟 ...
- 实验Matlab数值运算,MATLAB数值实验一(数据的插值运算及其应用完整版
<MATLAB数值实验一(数据的插值运算及其应用完整版>由会员分享,可在线阅读,更多相关<MATLAB数值实验一(数据的插值运算及其应用完整版(6页珍藏版)>请在人人文库网上搜 ...
- matlab bdir 排序,matlab-3次b样条（matlab - 3次b样条）.doc
matlab-3次b样条(matlab - 3次b样条) matlab-3次b样条(matlab - 3次b样条) http: / / / 3caddesign high speed, high pr ...
- matlab三点求法向量,matlab求法向量
|dT/ds| dT/ds 投影方向单位向量,垂直于 T 平面 T 和 N 的单位法向量,即曲率的平面 曲线的扭率: |dB/ds| 重力常数 力学中力的标准符号 弹簧的弹簧常数 ...... |dT ...
- B-spline三次B样条曲线方程
一.B-样条基函数 它有两条贝塞尔基函数所没有的特性, (1)定义域被节点细分(subdivided): (2) 基函数不是在整个区间非零.实际上,每个B样条基函数在附近一个子区间非零, 因此,B-样 ...
- 用gismo中B样条基函数替代自己写的基函数
文章目录 前言 一.需要添加[gismo库](https://blog.csdn.net/mw_1422102031/article/details/128966345?spm=1001.2014.3 ...
- 三次样条函数(cubic spline functions)的插值求解(python,数值积分)
第三十七篇 三次样条函数的插值求解 利用三次样条函数进行插值 到目前为止,描述的两种插值方法拉格朗日多项式和正向差分会形成高阶的多项式,但一般情况下,选择阶数的值等于小于数据点的数目更合适.除了考虑计 ...
最新文章
- .NET编码解码(HtmlEncode与HtmlEncode)
- html左边动右边不动,网页布局//上左不动,其他滑动
- python3 32位_Python 3.6.8软件安装教程
- 写出float x 与“零值”比较的if语句——一道面试题分析
- Nuget 启用数据库迁移的时候一定要把包含DbContext的项目设为启动项目
- code block怎样实现图形界面_微服务入门：Openresty实现API网关
- dj鲜生-36-商品应用-其它模型类的创建-完善goods应用的数据表
- html5移动端开发(rem和媒体查询@media)
- 转载 WebService 的CXF框架 WS方式Spring开发
- JavaScript学习（四十）—字面量创建对象图解
- html+css做的丝带标签
- AX 2009 父窗体参数记录传递
- 微软MVP申请“自我介绍”部分英文示范
- SQL必知必会 附录解读
- 如何备考系统集成项目管理工程师？
- 【写博客常用】Word文档中怎么插入分隔线
- 文本分类概述（nlp）
- DQL 数据查询语⾔
- Zbrush curve Tube ,Move topologica笔刷的使用
- VC环境中获取窗体标题栏的位置和高度