数值分析相关实验代码
MATLAB数值实验一(数据的插值运算及其应用完整版)
佛山科学技术学院实 验 报 告课程名称 数值分析 实验项目 插值法与数据拟合 专业班级 机械工程 姓 名 余红杰 学 号 10 指导教师 陈剑 成 绩 日 期 月 日一、实验目的1、学会Lagrange 插值、牛顿插值和三次样条插值等基本插值方法；2、讨论插值的Runge 现象3、学会Matlab 提供的插值函数的使用方法，会用这些函数解决实际问题。
二、实验原理1、拉格朗日插值多项式2、牛顿插值多项式3、三次样条插值 三、实验步骤1、用MATLAB 编写独立的拉格朗日插值多项式函数2、用MATLAB 编写独立的牛顿插值多项式函数3、用MATLAB 编写独立的三次样条函数（边界条件为第一、二种情形）4、已知函数在下列各点的值为：根据步骤1,2,3编好的程序，试分别用4次拉格朗日多项式4()L x 、牛顿插值多项式4()P x 以及三次样条函数()S x (自然边界条件)对数据进行插值，并用图给出 ｛(,),0.20.08,0,1,2,,10i i i x y x i i =+=｝，4()L x 、4()P x 和()S x 。
5、在区间[-1,1]上分别取10,20n =用两组等距节点对龙格函数21(),(11)125f x x x=-≤≤+作多项式插值，对不同n 值，分别画出插值函数及()f x 的图形。
6、下列数据点的插值可以得到平方根函数的近似，在区间[0,64]上作图。
（1）用这9个点作8次多项式插值8()L x 。
（2）用三次样条（第一边界条件）程序求()S x 。
7、对于给函数21()125f x x =+在区间[-1,1]上取10.2(0,1,,10)i x i i =-+=，试求3次曲线拟合，试画出拟合曲线并打印出方程，与第5题的结果比较。
四、实验过程与结果：1、Lagrange 插值多项式源代码：function ya=lag(x,y,xa) %x 所有已知插值点 %y 插值点对应函数值 %xa 所求点，自变量 %ya 所求点插值估计量 ya=0; mu=1; %初始化%循环方式求L 系数，并求和： for i = 1:length(y) for j = 1:length(x) if i ~= jmu = mu * (xa - x(j) ) / ( x(i) - x(j) ); else continue end endya = ya + y(i) * mu ; mu = 1; end2、Newton 源代码：function ya = newton(x,y,xa) %x 所有已知插值点 %y 插值点对应函数值 %xa 所求点，自变量 %ya 所求点插值估计量 %建立系数零矩阵D 及初始化：D = zeros(length(x)-1);ya = y(1);xi = 1;%求出矩阵D，该矩阵第一行为牛顿插值多项式系数：for i=1:(length(x)-1)D(i,1) = (y(i+1) -y(i))/(x(i+1) -x(i));endfor j=2:(length(x)-1)for i=1:(length(x)-j)D(i,j) = (D(i+1,j-1) - D(i,j-1)) / (x(i+j) - x(i)); endend%xi为单个多项式（x-x(1))(x-x(2))...的值for i=1:(length(x)-1)for j=1:ixi = xi*(xa - x(j));endya = ya + D(1,i)*xi;xi = 1;end3、三次样条插值多项式（1）（第一边界条件）源代码：function y=yt1(x0,y0,f_0,f_n,x) _____________(1)%第一类边界条件下三次样条插值；%xi 所求点；%yi 所求点函数值；%x 已知插值点；%y 已知插值点函数值；%f_0左端点一次导数值；%f_n右端点一次导数值；n = length(x0);z = length(y0);h = zeros(n-1,1);k=zeros(n-2,1);l=zeros(n-2,1);S=2*eye(n);for i=1:n-1h(i)= x0(i+1)-x0(i);endfor i=1:n-2k(i)= h(i+1)/(h(i+1)+h(i));l(i)= 1-k(i);end%对于第一种边界条件：k = [1;k]; _______________________(2)l = [l;1]; _______________________(3)%构建系数矩阵S：for i = 1:n-1S(i,i+1) = k(i);S(i+1,i) = l(i);end%建立均差表：F=zeros(n-1,2);for i = 1:n-1F(i,1) = (y0(i+1)-y0(i))/(x0(i+1)-x0(i));endD = zeros(n-2,1);for i = 1:n-2F(i,2) = (F(i+1,1)-F(i,1))/(x0(i+2)-x0(i));D(i,1) = 6 * F(i,2);end%构建函数D：d0 = 6*(F(1,2)-f_0)/h(1); ___________(4)dn = 6*(f_n-F(n-1,2))/h(n-1); ___________(5)D = [d0;D;dn]; ______________(6)m= S\D;%寻找x所在位置，并求出对应插值：for i = 1:length(x)for j = 1:n-1if (x(i)<=x0(j+1))&(x(i)>=x0(j))y(i) =( m(j)*(x0(j+1)-x(i))^3)/(6*h(j))+...(m(j+1)*(x(i)-x0(j))^3)/(6*h(j))+...(y0(j)-(m(j)*h(j)^2)/6)*(x0(j+1)-x(i))/h(j)+... (y0(j+1)-(m(j+1)*h(j)^2)/6)*(x(i)-x0(j))/h(j) ; break;else continue;endendend（2）（自然边界条件）源代码：仅仅需要对上面部分标注的位置做如下修改:__(1):function y=yt2(x0,y0,x)__(2):k=[0;k]__(3):l=[l;0]__(4)+(5):删除—(6):D=[0:D:0]4、——————————————PS:另建了一个f方程文件，后面有一题也有用到。
数值分析实验误差分析.doc
c[i]=temp;//保存首行信息
}
//消去l=aik/akk
k=j;
while(k<n-1){
i=j;
b[0] =a[k+1][j];//保留第一个系数防止后面破坏
for(;i<n+1;i++)
{
//b[1]=b[0]*c[i]/a[j][j];
a[k+1][i]=b[0]*c[i]/a[j][j]-a[k+1][i];
int n;
float fenzi[Max],fenmu[Max][Max];
void getdata()
{
cin>>n;
for(int i=0;i<n+1;i++)
{
cin>>xi[i]>>yi[i];
}
cin>>value_x;
}
void init()，
{
int i,j;
for(i=0;i<n+1;i++)
Gauss(a);
Gauss(b);
}
void Gauss(float a[n][n+1])
{
int max=0;
int i,j,k;
float b[n+1];//临时
float c[n+1];//
float x[n];//求解的x集合
float temp;
int couts=0;
// showarray(a);
通过这次实验我对拉格朗日插值和牛顿插值的原理的认识变得更加的深刻，明白了跟多编写此程序时要注意的问题。
数值分析课程实验报告-二分法和牛顿迭代法
用二分法和牛顿迭代法求方程的根
算法名称用二分法和牛顿迭代法求方程的根
学科专业机械工程
作者姓名xxxxxx
作者学号xxxBiblioteka xx作者班级xxxxxxxx
xx大学
二〇一五年十二月
《数值分析》课程实验报告
实验名称
用二分法和牛顿迭代法求方程的根
成绩
一、问题背景
在科学研究与工程计算中，常遇到方程（组）求根问题。若干个世纪以来，工程师和数学家花了大量时用于探索求解方程（组），研究各种各样的方程求解方法。对于方程f(x)=0，当f(x)为线性函数时，称f(x)=0为线性方程；当f(x）为非线性函数时，称式f(x)=0为非线性方程。对于线性方程（组）的求解，理论与数值求法的成果丰富；对于非线性方程的求解，由于f(x）的多样性，尚无一般的解析解法。当f(x)为非线性函数时，若f(x)=0无解析解，但如果对任意的精度要求，设计迭代方程，数值计算出方程的近似解，则可以认为求根的计算问题已经解决，至少能够满足实际要求。
fx=subs(ff,x,xk);
fa=subs(ff,x,a);
k=k+1;
iffx==0
y(k)=xk;
break;
elseiffa*fx<0
b=xk;
else
a=xk;
end
y(k)=xk;
end
plot(y,'.-');
gridon
（2）牛顿迭代法程序：
functionx=newton(xx,n)
对二分法和牛顿迭代法的观察和分析我们可以知道，二分法的优点是方法比较简单，编程比较容易，只是二分法只能用于求方程的近似根，不能用于求方程的复根，且收敛速度慢。而牛顿迭代法的收敛速度明显大于二分法的速度。
数值计算基础实验报告(3篇)
第1篇一、实验目的1. 理解数值计算的基本概念和常用算法；2. 掌握Python编程语言进行数值计算的基本操作；3. 熟悉科学计算库NumPy和SciPy的使用；4. 分析算法的数值稳定性和误差分析。
二、实验内容1. 实验环境操作系统：Windows 10编程语言：Python 3.8科学计算库：NumPy 1.19.2，SciPy 1.5.02. 实验步骤（1）Python编程基础1）变量与数据类型2）运算符与表达式3）控制流4）函数与模块（2）NumPy库1）数组的创建与操作2）数组运算3）矩阵运算（3）SciPy库1）求解线性方程组2）插值与拟合3）数值积分（4）误差分析1）舍入误差2）截断误差3）数值稳定性三、实验结果与分析1. 实验一：Python编程基础（1）变量与数据类型通过实验，掌握了Python中变量与数据类型的定义方法，包括整数、浮点数、字符串、列表、元组、字典和集合等。
（2）运算符与表达式实验验证了Python中的算术运算、关系运算、逻辑运算等运算符，并学习了如何使用表达式进行计算。
（3）控制流实验学习了if-else、for、while等控制流语句，掌握了条件判断、循环控制等编程技巧。
（4）函数与模块实验介绍了Python中函数的定义、调用、参数传递和返回值，并学习了如何使用模块进行代码复用。
2. 实验二：NumPy库（1）数组的创建与操作通过实验，掌握了NumPy数组的基本操作，包括创建数组、索引、切片、排序等。
（2）数组运算实验验证了NumPy数组在数学运算方面的优势，包括加、减、乘、除、幂运算等。
（3）矩阵运算实验学习了NumPy中矩阵的创建、操作和运算，包括矩阵乘法、求逆、行列式等。
3. 实验三：SciPy库（1）求解线性方程组实验使用了SciPy库中的线性代数模块，通过高斯消元法、LU分解等方法求解线性方程组。
（2）插值与拟合实验使用了SciPy库中的插值和拟合模块，实现了对数据的插值和拟合，并分析了拟合效果。
学科专业代码081200
秋季
S7700015Q
高级运筹学
2
36
秋季
S7700038Q
小波分析
2
36
秋季
S7700050C
网络科学导论
1
18
春季
S7700051Q
复杂网络建模
2
36
秋季
专业基础课
S7101001Q
人工智能
2
36
秋季
S7101011C
数据库系统原理
2
36
春季
S7101002Q
算法设计与分析
2
28/8
秋季
S7101003Q
S7700005QN
历史文明与国际关系
1
18
秋季
S7700047Q
幸福与公正
1
18
秋季
S7700048Q
文艺与审美鉴赏
1
18
秋季
必修环节
学术活动
1～2学分(学术活动参加5次～9次，记为1学分；参加10次或10次以上记为2学分)
学科专业代码：081200
学科专业名称：计算机科学与技术
类型：学术研究型
研究方向：
1.计算机应用技术2.计算机软件与理论
3.模式识别与智能系统4.信息安全
课程设置：
课程性质
类别
课程编号
课程名称
学分
学时
课内/实验
学期
备注
学
位
课
通识课
S7700001Q
中国特色社会主义理论与实践研究
2
36
秋季
S7700002Q
计算分子生物学
2
36
秋季
0～4学分为跨一级学科专业课
数值分析实验2014
数值分析实验(2014,9,16~10,28)信计1201班，人数34人数学系机房数值分析计算实习报告册专业__________________学号_______________姓名_______________2014~2015年第一学期实验一数值计算的工具Matlab1. 解释下MATLABS序的输出结果程序：t=0.1n=1:10e=n/10-n*te 的结果：0 0 -5.5511e-017 0 0-1.1102e-016 -1.1102e-016 0 0 02. 下面MATLABS序的的功能是什么？程序：x=1;while 1+x>1,x=x/2,pause(0.02),e nd用迭代法求出x=x/2,的最小值x=1;while x+x>x,x=2*x,pause(0.02),e nd用迭代法求出x=2*x,的值，使得2x>Xx=1;while x+x>x,x=x/2,pause(0.02),e nd用迭代法求出x=x/2,的最小值，使得2x>X3. 考虑下面二次代数方程的求解问题2ax bx c = 0公式x=电上4ac是熟知的，与之等价地有_____________________________ ,对于2a-b ■ b -4aca =1,b =100000000,c =1，应当如何选择算法。
b ~4ac计算，因为b与b2— 4ac相近，两个相加减不宜应该用2a u做分母3 5 74. 函数sin(x)有幂级数展开sin x = x - x - - ■■3! 5! 7!利用幕级数计算sinx的MATLAB程序为fun cti on s=powers in(x)s=0;t=x;n=1;while s+t~=s;s=s+t ；t=-x A2/ ((n+1)*(n+2) ) *t ;n=n+2 ；endt仁cputime;pause(10);t2=cputime;t0=t2-t1(a) 解释上述程序的终止准则。
数值分析实验报告
%消元过程
fori=k+1:n
m=A(i,k)/A(k,k);
forj=k+1:n
A(i,j)=A(i,j)-m*A(k,j);
end
b(i)=b(i)-m*b(k);
end
det=det*A(k,k);
end
det=det*A(n,n);
%回代过程
ifabs(A(n,n))<1e-10
flag='failure';return;
*x=(x0,x1….,xn),插值节点
*y=(y0,y1,…,yn);被插函数f(x)在插值节点处的函数值
*t求插值函数Pn(x)在t处的函数值
*返回值 插值函数Pn(x)在t处的函数值
*/
procedureNewton
forj=0to n
d1jyj;
endfor
forj=1to n
fori=j to n
[n,m]=size(A);nb=length(b)
%当方程组行与列的维数不相等时，停止计算，并输出出错信息
ifn~=m
error('The row and columns of matrix A must beepual!');
return;
end
%当方程组与右端项的维数不匹配时，停止计算，并输出错误信息
clear
fprintf('gauss-seidel迭代法')
x1_(1)=0;
x2_(1)=0;
x3_(1)=0;
fori=1:9
x1_(i+1)=7.2+0.1*x2_(i)+0.2*x3_(i);
数值分析方法计算定积分
数值分析⽅法计算定积分⽤C语⾔实现⼏种常⽤的数值分析⽅法计算定积分，代码如下：1 #include <stdio.h>2 #include <string.h>3 #include <stdlib.h>4 #include <math.h>56#define EPS 1.0E-14 //计算精度7#define DIV 1000 //分割区间，值越⼤，精确值越⾼8#define ITERATE 20 //⼆分区间迭代次数，区间分割为2^n，初始化应该⼩⼀点，否则会溢出910#define RECTANGLE 0 //矩形近似11#define TRAPEZOID 1 //梯形近似12#define TRAPEZOID_FORMULA 2 //递推梯形公式13#define SIMPSON_FORMULA 3 //⾟普森公式14#define BOOL_FORMULA 4 //布尔公式1516double function1(double);17double function2(double);18double function3(double);19void Integral(int, double f(double), double, double); //矩形, 梯形逼近求定积分公式20void Trapezoid_Formula(double f(double), double a, double b); //递推梯形公式21void Simpson_Formula(double f(double), double a, double b); //⾟普森公式22void Bool_Formula(double f(double), double a, double b); //布尔公式23void Formula(int, double f(double), double, double);2425double function1(double x)26 {27double y;28 y = cos(x);29return y;30 }3132double function2(double x)33 {34double y;35 y=1/x;36// y = 2+sin(2 * sqrt(x));37return y;38 }3940double function3(double x)41 {42double y;43 y = exp(x);44return y;45 }4647int main()48 {49 printf("y=cos(x) , x=[0, 1]\n");50 printf("Rectangle : "); Integral(RECTANGLE, function1, 0, 1);51 printf("Trapezoid : "); Integral(TRAPEZOID, function1, 0, 1);52 Formula(TRAPEZOID_FORMULA, function1, 0, 1);53 Formula(SIMPSON_FORMULA, function1, 0, 1);54 Formula(BOOL_FORMULA, function1, 0, 1);55 printf("\n");5657 printf("y=1/x , x=[1, 5]\n");58 printf("Rectangle : "); Integral(RECTANGLE, function2, 1, 5);59 printf("Trapezoid : "); Integral(TRAPEZOID, function2, 1, 5);60 Formula(TRAPEZOID_FORMULA, function2, 1, 5);61 Formula(SIMPSON_FORMULA, function2, 1, 5);62 Formula(BOOL_FORMULA, function2, 1, 5);63 printf("\n");6465 printf("y=exp(x) , x=[0, 1]\n");66 printf("Rectangle : "); Integral(RECTANGLE, function3, 0, 1);67 printf("Trapezoid : "); Integral(TRAPEZOID, function3, 0, 1);68 Formula(TRAPEZOID_FORMULA, function3, 0, 1);69 Formula(SIMPSON_FORMULA, function3, 0, 1);70 Formula(BOOL_FORMULA, function3, 0, 1);7172return0;73 }74//积分：分割-近似-求和-取极限75//利⽤黎曼和求解定积分76void Integral(int method, double f(double),double a,double b) //f(double)，f(x), double a，float b，区间两点77 {79double x;80double sum = 0; //矩形总⾯积81double h; //矩形宽度82double t = 0; //单个矩形⾯积8384 h = (b-a) / DIV;8586for(i=1; i <= DIV; i++)87 {88 x = a + h * i;89switch(method)90 {91case RECTANGLE :92 t = f(x) * h;93break;94case TRAPEZOID :95 t = (f(x) + f(x - h)) * h / 2;96break;97default:98break;99 }100 sum = sum + t; //各个矩形⾯积相加101 }102103 printf("the value of function is %.10f\n", sum);104 }105106//递推梯形公式107void Trapezoid_Formula(double f(double), double a, double b)//递推梯形公式108 {109 unsigned int i, j, M, J=1;110double h;111double F_2k_1 = 0;112double T[32];113double res = 0;114115 T[0] = (f(a) + f(b)) * (b-a)/2;116for(i=0; i<ITERATE; i++)117 {118 F_2k_1 = 0;119 J = 1;120121for(j=0; j<=i; j++) //区间划分122 J <<= 1; //2^J123 h = (b - a) /J; //步长124//M = J/2; //2M表⽰将积分区域划分的块数125for(j=1; j <= J; j += 2) //126 {127 F_2k_1 += f(a + j*h); //f(2k-1)累加和128 }129 T[i+1] = T[i]/2 + h * F_2k_1; //递推公式130 res = T[i+1];131132if(fabs(T[i+1] - T[i]) < EPS)133if(i < 3) continue;134else break;135 }136137 printf("Trapezoid_Formula : the value of function is %.10f\n", res);138//return res;139 }140//⾟普森公式141void Simpson_Formula(double f(double), double a, double b) //⾟普森公式142 {143 unsigned int i, j, M, J=1;144double h;145double F_2k_1 = 0;146double T[32], S[32];147double res_T = 0, res_S = 0;148149 T[0] = (f(a) + f(b)) * (b-a)/2;150for(i=0; i<ITERATE; i++)151 {152 F_2k_1 = 0;153 J = 1;154155for(j=0; j<=i; j++) //区间划分156 J <<= 1; //2^J157 h = (b - a) /J; //步长158//M = J/2; //2M表⽰将积分区域划分的块数159for(j=1; j <= J; j += 2) //160 {161 F_2k_1 += f(a + j*h); //f(2k-1)累加和163 T[i+1] = T[i]/2 + h * F_2k_1; //递推梯形公式164 res_T = T[i+1];165166if(fabs(T[i+1] - T[i]) < EPS)167if(i < 3) continue;168else break;169 }170171for(i=1; i<ITERATE; i++)172 {173 S[i] = (4 * T[i] - T[i-1]) / 3; //递推⾟普森公式174 res_S = S[i];175if(fabs(S[i] - S[i-1]) < EPS)176if(i < 3) continue;177else break;178 }179180 printf("Simpson_Formula : the value of function is %.10f\n", res_S);181//return res_S;182 }183184//布尔公式185void Bool_Formula(double f(double), double a, double b) //布尔公式186 {187 unsigned int i, j, M, J=1;188double h;189double F_2k_1 = 0;190double T[32], S[32], B[32];191double res_T = 0, res_S = 0, res_B;192193 T[0] = (f(a) + f(b)) * (b-a)/2;194for(i=0; i<ITERATE; i++)195 {196 F_2k_1 = 0;197 J = 1;198199for(j=0; j<=i; j++) //区间划分200 J <<= 1; //2^J201 h = (b - a) /J; //步长202//M = J/2; //2M表⽰将积分区域划分的块数203for(j=1; j <= J; j += 2) //204 {205 F_2k_1 += f(a + j*h); //f(2k-1)累加和206 }207 T[i+1] = T[i]/2 + h * F_2k_1; //递推公式208 res_T = T[i+1];209210if(fabs(T[i+1] - T[i]) < EPS)211if(i < 3) continue;212else break;213 }214215for(i=1; i<ITERATE; i++)216 {217 S[i] = (4 * T[i] - T[i-1]) / 3; //递推⾟普森公式218 res_S = S[i];219if(fabs(S[i] - S[i-1]) < EPS)220if(i < 3) continue;221else break;222 }223224for(i=1; i<ITERATE; i++)225 {226 B[i] = (16 * S[i] - S[i-1]) / 15; //递推布尔公式227 res_B = B[i];228if(fabs(B[i] - B[i-1]) < EPS)229if(i < 3) continue;230else break;231 }232233 printf("Bool_Formula : the value of function is %.10f\n", res_B);234//return res_B;235 }236237//采⽤⼆分区间的⽅法迭代，实际运⾏速度很慢238void Formula(int method, double f(double), double a, double b)//递推梯形公式239 {240 unsigned int i, j, M, J=1;241double h;242double F_2k_1 = 0;243double T[32], S[32], B[32];244double res_B = 0, res_S = 0, res_T = 0;247for(i=0; i<ITERATE; i++)248 {249 F_2k_1 = 0;250 J = 1;251252for(j=0; j<=i; j++) //区间划分253 J <<= 1; //2^J254 h = (b - a) /J; //步长255//M = J/2; //2M表⽰将积分区域划分的块数256for(j=1; j <= J; j += 2) //257 {258 F_2k_1 += f(a + j*h); //f(2k-1)累加和259 }260 T[i+1] = T[i]/2 + h * F_2k_1; //递推公式261 res_T = T[i+1];262263if(fabs(T[i+1] - T[i]) < EPS)264if(i < 3) continue;265else break;266 }267268switch(method)269 {270default :271case TRAPEZOID_FORMULA :272 printf("Trapezoid_Formula : the value of function is %.10f\n", res_T); 273//return res_T;274break;275case SIMPSON_FORMULA :276for(i=1; i<ITERATE; i++)277 {278 S[i] = (4 * T[i] - T[i-1]) / 3;279 res_S = S[i];280if(fabs(S[i] - S[i-1]) < EPS)281if(i < 3) continue;282else break;283 }284 printf("Simpson_Formula : the value of function is %.10f\n", res_S); 285//return res_S;286break;287case BOOL_FORMULA :288for(i=1; i<ITERATE; i++)289 {290 S[i] = (4 * T[i] - T[i-1]) / 3;291 res_S = S[i];292if(fabs(S[i] - S[i-1]) < EPS)293if(i < 3) continue;294else break;295 }296for(i=1; i<ITERATE; i++)297 {298 B[i] = (16 * S[i] - S[i-1]) / 15;299 res_B = B[i];300if(fabs(B[i] - B[i-1]) < EPS)301if(i < 3) continue;302else break;303 }304305 printf("Bool_Formula : the value of function is %.10f\n", res_B); 306//return res_B;307break;308 }309310 }测试结果：。
三次样条插值matlab代码实现
三次样条插值matlab代码实现
三次样条插值是一种常用的数值分析方法，用于在给定的数据点上拟合出一个光滑的曲线。
在Matlab中，可以使用内置的spline函数来实现三次样条插值。
以下是一个简单的示例代码：
matlab.
% 创建一些示例数据点。
x = 1:5;
y = [3 6 5 8 9];
% 使用spline函数进行三次样条插值。
xx = 1:0.1:5;
yy = spline(x, y, xx);
% 绘制原始数据点和插值结果。
plot(x, y, 'o', xx, yy, '-');
legend('原始数据', '插值结果');
在这个示例中，我们首先创建了一些示例数据点x和y。
然后使用spline函数对这些数据点进行三次样条插值，得到了插值结果xx和yy。
最后，我们使用plot函数将原始数据点和插值结果进行了可视化展示。
需要注意的是，样条插值是一种较为复杂的数值计算方法，需要对输入数据进行适当的处理和理解。
在实际应用中，可能需要根据具体情况对插值方法进行调整和优化，以获得更好的结果。
希望这个简单的示例能够帮助你理解如何在Matlab中实现三次样条插值。
如果你有更多的问题或者需要进一步的解释，请随时告诉我。
数值分析实验复合辛普森公式报告
实验题目：用复合辛普森公式求方程的积学生姓名：***专业：信息与计算科学学号：**********完成日期：2011/11/20西安科技大学计算机科学与技术学院实验题目：利用复合辛普森公式求方程⎰-=109/4ln xdx x 的积学生姓名： 何弯弯 学号： 0908060222 完成日期：2011/10/22 1 实验目的 1、了解复合辛普森公式的方法原理；2、利用复合辛普森公式求方程⎰-=109/4ln xdx x 的积2.1 算法原理原理：将区间[a,b]等分成N 个子区间[x(k),x(k+1)]（k=0,1,…..,N -1），h=（b-a ）/N ，在每个子区间[x(k),x(k+1)]上使用辛普森公式，可求得结果。
辛普森公式:])(4)(2)()([6/)]()()([6/1115.01015.0∑∑∑=-=+-=+++++=++=n k n k k k n k k k k n x f x f b f a f h x f x f x f h S其中5.0+x x 是[1,+k k x x ]得中点，即2/5.0h x x k k +=+2.2 算法步骤步骤一：输入a,b,h步骤二：代入辛普森公式求的结果2.3 程序流程图3 实验结果分析4实验心得体会通过本次实验我熟悉了用复化辛普森公式求数值积分的全过程，并更加熟悉的掌握了复化辛普生公式以及原理。
使我对数值分析这门课有了进一步的理解，为以后的学习打下了良好的基础。
参考文献[1]龚尚福，贾彭涛，靳玉萍.C/C++语言程序设计.徐州：中国矿业大学出版社，2006[2]李庆阳，王能超，易大义.数值分析（第五版）.北京：清华大学出版社，2008附录（源代码）复合辛普森算法代码：#include<iostream.h>#include<math.h>double get(double x){if(x==0) x=0.000000000001;return sqrt(x)*log(x);}double g(double a,double h,double b){int i,n;double sum1,sum2=0,f;n=(int)((b-a)/h);sum1=get(a+h/2);for(i=1;i<=n-1;i++) //求f(x(k))的和{ sum1+=get(a+i*h+h/2);sum2+=get(a+i*h);}f=4*sum1+2*sum2;return f;}void main(){ double a=0,b=1,h;//给定的区间，以及步长int n;while(true){cout<<"请输入要给的步长h(h>0):";cin>>h;if(h==-1)break;double I1,f1,f2,f3;n=(int)((b-a)/h);f1=get(a);cout<<f1<<endl;f2=g(a,h,b);cout<<f2<<endl;f3=get(b);cout<<f3<<endl;I1=h/6*(f1+f2+f3);//复化辛普森公式cout<<"给定步长为h="<<h<<endl;cout<<"复化辛普森公式计算结果:I="<<I1<<endl;cout<<"所分成n个区间，n="<<n<<endl;}}。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
连续复合梯形和辛普森公式：
function y=hanshu(x)
y=(2*exp(-x))/sqrt(pi);
………………………………
x=0:0.01:1;
h=0.01;
futix=0;
for n=1:100
y1=hanshu(x(n));
y2=hanshu(x(n+1));
futix=futix+(y1+y2)*h/2;
end
futix=vpa(futix,6)
fusim=0;
for n=1:100
y1=hanshu(x(n));
y2=hanshu(x(n+1));
y3=hanshu((x(n)+x(n+1))/2);
fusim=fusim+(y1+y2+4*y3)*h/6;
end
fusim=vpa(fusim,6)
离散复合梯形和辛普森公式：
x=[0,1,4,9,16,25,36,49,64];
y=[0,1,2,3,4,5,6,7,8];
plot(x,y,'ro')
x1=0:0.1:64;
y1=interp1(x,y,x1,'spline');
hold on
plot(x1,y1,'b.')
hold off
ftix=0;
h=0.1;
for n=1:length(x1)-1
ftix=ftix+(y1(n)+y1(n+1))*h/2;
end
ftix=vpa(ftix,10)
fsim=0;
h=h*2;
for n=1:2:length(x1)-2
fsim=fsim+(y1(n)+y1(n+2)+4*y1(n+1))*h/6;
end
fsim=vpa(fsim,10)
改进的欧拉方法：
functiondy=hanshu(x,y)
dy=1/x^2-y/x;
………………………………
x=1:0.05:2;
y(1)=1;
h=0.05;
for n=1:length(x)-1
y1=y(n)+h*hanshu(x(n),y(n));
y(n+1)=y(n)+h*(hanshu(x(n),y(n)))+hanshu((x(n)+h),y1)/2;
end
y
龙格库塔方法：
functiondy=h(t,y)
dy=zeros(3,1);
dy(1)=-y(1);
dy(2)=y(1)-10*y(2);
dy(3)=-10*y(2);
…………………………
tspan=0:0.01:1;
[T,Y] = ode45(@h,tspan,[1 1 1]);
plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')
LU分解：
A=[1,2,3;2,5,2;3,1,5];
b=[14;18;20];
ifdet(A)==0
error('');
else
[m,n]=size(A);
L=zeros(m,m);
L=L+diag(ones(m,1));
U=zeros(m,m);
U(1,:)=A(1,:);
L(2:end,1)=A(2:end,1)/U(1,1);
for r=2:m
U(r,r)=A(r,r)-L(r,1:r-1)*U(1:r-1,r);
for i=r+1:m
U(r,i)=A(r,i)-L(r,1:r-1)*U(1:r-1,i);
L(i,r)=(A(i,r)-L(i,1:r-1)*U(1:r-1,r))/U(r,r);
end
end
end
L
U
y(1)=b(1);
for i=2:n
y(i)=b(i)-L(i,1:i-1)*y(1:i-1)';
end
x(n)=y(n)/U(n,n);
for i=n-1:-1:1
x(i)=(y(i)-U(i,i+1:n)*x(i+1:n)')/U(i,i);
end
x