龙格库塔方法matlab实现
龙格库塔方法matlab实现~ function ff=rk(yy,x0,y0,h,a,b)%yy为y的导函数,x0,y0,为初值,h为步长,a,b为区间 c=(b-a)/h+1;i1=1; %c为迭代步数;i1为迭代步数累加值 y=y0;z=zeros(c,6); %z生成c行,5列的零矩阵存放结果; %每行存放c次迭代结果,每列分别存放k1~k4及y的结果 for x=a:h:b if i1<=c k1=feval(yy,x,y); k2=feval(yy,x+h/2,y+(h*k1)/2); k3=feval(yy,x+h/2,y+(h*k2)/2); k4=feval(yy,x+h,y+h*k3); y=y+(h/6)*(k1+2*k2+2*k3+k4); z(i1,1)=x;z(i1,2)=k1;z(i1,3)=k2;z(i1,4)=k3;z(i1,5)=k4;z(i1,6)=y; i1=i1+1; end end fprintf(‘结果矩阵,第一列为x(n),第二列~第五列为k1~k4,第六列为y(n+1)的结果') z %在命令框输入下列语句 %yy=inline('x+y'); %>> rk(yy,0,1,0.2,0,1) %将得到结果 %结果矩阵,第一列为x(n),第二列~第五列为k1~k4第六列为y(n+1)的结果 %z = % 0 1.0000 1.2000 1.2200 1.4440 1.2428 % 0.2000 1.4428 1.6871 1.7115 1.9851 1.5836 % 0.4000 1.9836 2.2820 2.3118 2.6460 2.0442 % 0.6000 2.6442 3.0086 3.0451 3.4532 2.6510 % 0.8000 3.4510 3.8961 3.9407 4.4392 3.4365 % 1.0000 4.4365 4.9802 5.0345 5.6434 4.4401
拉格朗日插值法C语言的实现
实验 一 .拉格朗日插值法C 语言的实现 1.实验目的: 进一步熟悉拉格朗日插值法。 掌握编程语言字符处理程序的设计和调试技术。 2.实验要求: 已知:某些点的坐标以及点数。 输入:条件点数以及这些点的坐标 。 输出:根据给定的点求出其对应的拉格朗日插值多项式的值 。 3.程序流程: (1)输入已知点的个数; (2)分别输入已知点的X 坐标; (3)分别输入已知点的Y 坐标; (4)通过调用函数lagrange 函数,来求某点所对应的函数值。 拉格朗日插值多项式如下: 0L ()()0,1,n n j k k j j k x y l x y j n ====∑…… 其中00()()0,1,,()k k x x l x k n x x -= =-k-1k+1n k k-1k k+1k n ……(x-x )(x-x )?…(x-x )…………(x -x )(x -x )?…(x -x ) 程序流程图:
↓ 程序如下: #include #include <> #include <> float lagrange(float *x,float *y,float xx,int n) /*拉格朗日插值算法*/ { int i,j; float *a,yy=; /*a作为临时变量,记录拉格朗日插值多项式*/ a=(float *)malloc(n*sizeof(float)); for(i=0;i<=n-1;i++) { a[i]=y[i]; for(j=0;j<=n-1;j++) if(j!=i) a[i]*=(xx-x[j])/(x[i]-x[j]); yy+=a[i]; } free(a); return yy; } int main() { int i; int n; float x[20],y[20],xx,yy; printf("Input n:"); scanf("%d",&n); if(n>=20) { printf("Error!The value of n must in (0,20)."); getch();return 1; } if(n<=0) { printf("Error! The value of n must in (0,20)."); getch(); return 1; } for(i=0;i<=n-1;i++) { printf("x[%d]:",i); scanf("%f",&x[i]); } printf("\n");
数值分析拉格朗日插值法上机实验报告
课题一:拉格朗日插值法 1.实验目的 1.学习和掌握拉格朗日插值多项式。 2.运用拉格朗日插值多项式进行计算。 2.实验过程 作出插值点(1.00,0.00),(-1.00,-3.00),(2.00,4.00)二、算法步骤 已知:某些点的坐标以及点数。 输入:条件点数以及这些点的坐标。 输出:根据给定的点求出其对应的拉格朗日插值多项式的值。 3.程序流程: (1)输入已知点的个数; (2)分别输入已知点的X坐标; (3)分别输入已知点的Y坐标; 程序如下: #include #include #include float lagrange(float *x,float *y,float xx,int n) /*拉格朗日
插值算法*/ { int i,j; float *a,yy=0.0; /*a作为临时变量,记录拉格朗日插值多项*/ a=(float*)malloc(n*sizeof(float)); for(i=0;i<=n-1;i++) { a[i]=y[i]; for(j=0;j<=n-1;j++) if(j!=i) a[i]*=(xx-x[j])/(x[i]-x[j]); yy+=a[i]; } free(a); return yy; } int main() { int i; int n; float x[20],y[20],xx,yy; printf("Input n:");
scanf("%d",&n); if(n<=0) { printf("Error! The value of n must in (0,20)."); getch();return 1; } for(i=0;i<=n-1;i++) { printf("x[%d]:",i); scanf("%f",&x[i]); } printf("\n"); for(i=0;i<=n-1;i++) { printf("y[%d]:",i);scanf("%f",&y[i]); } printf("\n"); printf("Input xx:"); scanf("%f",&xx); yy=lagrange(x,y,xx,n); printf("x=%f,y=%f\n",xx,yy); getch(); } 举例如下:已知当x=1,-1,2时f(x)=0,-3,4,求f(1.5)的值。
matlab实现Lagrange多项式插值观察龙格现象
Matlab进行Lagrange多项式插值 拉格朗日插值法对函数y=1./(1+25*x.^2)在区间[-1,1]进行5次、10次、15次插值观察龙格现象 主程序 1.拉格朗日 function [c,l]=lagran(x,y) %c为多项式函数输出的系数 %l为矩阵的系数多项式 %x为横坐标上的坐标向量 %y为纵坐标上的坐标向量 w=length(x); n=w-1; l=zeros(w,w); for k=1:n+1 v=1; for j=1:n+1 if k~=j v=conv(v,poly(x(j)))/(x(k)-x(j)) %对多项式做卷积运算 end end l(k,:)=v; end c=y*l; 2.在matlab窗口中输入: x=linspace(-1,1,6);y=1./(1+25*x.^2); lagran(x,y) 回车可得结果: ans = -0.0000 1.2019 -0.0000 -1.7308 -0.0000 0.5673 在matlab窗口中输入: x=linspace(-1,1,11);y=1./(1+25*x.^2); lagran(x,y) 回车可得结果: ans = -220.9417 0.0000 494.9095 -0.0000 -381.4338 -0.0000 123.3597 0.0000 -16.8552 0.0000 1.0000 在matlab窗口中输入: x=linspace(-1,1,16);y=1./(1+25*x.^2); lagran(x,y) 回车可得结果: ans =
1.0e+003 * Columns 1 through 14 0.0000 -1.5189 -0.0000 4.6511 -0.0000 -5.5700 0.0000 3.3477 0.0000 -1.0830 -0.0000 0.1901 -0.0000 -0.0180 Columns 15 through 16 0.0000 0.0010 3.由以上结果可定义一下函数: function y=f1(x) y=1./(1+25*x.^2); function y=f2(x) y=1.2019*x.^4 -1.7308*x.^2+0.5673; function y=f3(x) y=-220.9417*x.^10+494.9095*x.^8-381.4338*x.^6+123.3597*x.^4-16.8552*x.^2+1; function y=f4(x) y=1*10^3*(-1.5189*x.^14+4.6511*x.^12-5.5700*x.^10+3.3477*x.^8-1.0830*x.^6+0.1901*x.^4-0.0180*x.^2+0.0010) 4. 在matlab窗口中输入: s1=@f1;s2=@f2;s3=@f3;s4=@f4;fplot(s1,[-1 1],'r');hold on;fplot(s2,[-1 1],'k');hold on;fplot(s3,[-1 1],'g');hold on;fplot(s4,[-1 1],'b');xlabel('input');ylabel('output');title('龙格现象');legend('s1=f(x)','s2=L5(x)','s3=L10(x)','s4=L15(X)');grid on 可以得到下图:
Matlab中龙格-库塔(Runge-Kutta)方法原理及实现
函数功能编辑本段回目录 ode是专门用于解微分方程的功能函数,他有ode23,ode45,ode23s等等,采用的是Runge-Kutta算法。ode45表示采用四阶,五阶runge-kutta单步算法,截断误差为(Δx)3。解决的是Nonstiff(非刚性)的常微分方程.是解决数值解问题的首选方法,若长时间没结果,应该就是刚性的,换用ode23来解. 使用方法编辑本段回目录 [T,Y] = ode45(odefun,tspan,y0) odefun 是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名 tspan 是区间[t0 tf] 或者一系列散点[t0,t1,...,tf] y0 是初始值向量 T 返回列向量的时间点 Y 返回对应T的求解列向量 [T,Y] = ode45(odefun,tspan,y0,options) options 是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等 [T,Y,TE,YE,IE] =ode45(odefun,tspan,y0,options) 在设置了事件参数后的对应输出 TE 事件发生时间 YE 事件解决时间 IE 事件消失时间 sol =ode45(odefun,[t0 tf],y0...) sol 结构体输出结果 应用举例编辑本段回目录 1 求解一阶常微分方程
程序: 一阶常微分方程 odefun=@(t,y) (y+3*t)/t^2; %定义函数 tspan=[1 4]; %求解区间 y0=-2; %初值 [t,y]=ode45(odefun,tspan,y0); plot(t,y) %作图 title('t^2y''=y+3t,y(1)=-2,1插值法实验报告
实验二插值法 1、实验目的: 1、掌握直接利用拉格郎日插值多项式计算函数在已知点的函数值;观察拉格郎日插值的龙格现象。 2、了解Hermite插值法、三次样条插值法原理,结合计算公式,确定函数值。 2、实验要求: 1)认真分析题目的条件和要求,复习相关的理论知识,选择适当的解决方案和算法; 2)编写上机实验程序,作好上机前的准备工作; 3)上机调试程序,并试算各种方案,记录计算的结果(包括必要的中间结果); 4)分析和解释计算结果; 5)按照要求书写实验报告; 3、实验内容: 1) 用拉格郎日插值公式确定函数值;对函数f(x)进行拉格郎日插值,并对f(x)与插值多项式的曲线作比较。 已知函数表:(0.56160,0.82741)、(0.56280,0.82659)、(0.56401,0.82577)、(0.56521,0.82495)用三次拉格朗日插值多项式求x=0.5635时函数近似值。 2) 求满足插值条件的插值多项式及余项 1) 4、题目:插值法 5、原理: 拉格郎日插值原理: n次拉格朗日插值多项式为:L n (x)=y l (x)+y 1 l 1 (x)+y 2 l 2 (x)+…+y n l n (x)
n=1时,称为线性插值, L 1(x)=y (x-x 1 )/(x -x 1 )+y 1 (x-x )/(x 1 -x )=y +(y 1 -x )(x-x )/(x 1 -x ) n=2时,称为二次插值或抛物线插值, L 2(x)=y (x-x 1 )(x-x 2 )/(x -x 1 )/(x -x 2 )+y 1 (x-x )(x-x 2 )/(x 1 -x )/(x 1 -x 2 )+y 2 (x -x 0)(x-x 1 )/(x 2 -x )/(x 2 -x 1 ) n=i时, Li= (X-X0)……(X-X i-1)(x-x i+1) ……(x-x n) (X-X0)……(X-X i-1)(x-x i+1) ……(x-x n) 6、设计思想: 拉格朗日插值法是根据n + 1个点x0, x1, ... x n(x0 < x1 < ... x n)的函数值f (x0), f (x1) , ... , f (x n)推出n次多項式p(x),然后n次多項式p (x)求出任意的点x对应的函数值f (x)的算法。 7、对应程序: 1 ) 三次拉格朗日插值多项式求x=0.5635时函数近似值 #include"stdio.h" #define n 5 void main() { int i,j; float x[n],y[n]; float x1; float a=1; float b=1; float lx=0; printf("\n请输入想要求解的X:\n x="); scanf("%f",&x1); printf("请输入所有点的横纵坐标:\n"); for(i=1;i实验二 拉格朗日插值龙格现象
汕 头 大 学 实 验 报 告 学院: 工学院系: 计算机系专业: 计算机科学与技术年级:2010 姓名: 林金正学号:2010101032完成实验时间: 5月24日 一.实验名称:拉格朗日插值的龙格现象 二.实验目的: 通过matlab 处理,观察拉格朗日插值的龙格现象. 三.实验内容: (1)学习matlab 的使用 (2)以实验的方式,理解高阶插值的病态性,观察拉格朗日插值的龙格现象。 四.实验时间、地点,设备: 实验时间:5月24日 实验地点:宿舍 实验设备:笔记本电脑 五,实验任务 在区间[-5,5]上取节点数n=11,等距离h=1的节点为插值点,对于函数2 5()1f x x =+进行拉格朗日插值,把f(x)与插值多项式的曲线花在同一张图上。 六.实验过程 拉格朗日插值函数定义: 对某个多项式函数,已知有给定的k + 1个取值点: 其中对应著自变数的位置,而对应著函数在这个位置的取值。 假设任意两个不同的xj 都互不相同,那麼应用拉格朗日插值公式所得到的拉格朗日插值多项式为: 其中每个为拉格朗日基本多项式(或称插值基函数),其表达式为: [3] 拉格朗日基本多项式 的特点是在上取值为1,在其它的点上取值为0。
1.使用matlab,新建function.m 文件,使用老师所给代码,构建拉格朗日函数. %lagrange.m function y=lagrange(x0,y0,x) n=length(x0); m=length(x); fori=1:m z=x(i);s=0; for k=1:n L=1; for j=1:n if j~=k L=L*(z-x0(j))/(x0(k)-x0(j)); end end s=s+L*y0(k); end y(i)=s; end y; 程序解释: (x0,y0):已知点坐标 x:所求点的横坐标, y:由(x0,y0)所产生的插值函数,以x 为参数,所的到的值 2.再一次新建function.m 文件. 构建自定义函数:25()1f x x = + %f.m function y = f(x) y = 5/(1+x*x); end 3.在脚本窗口中输入: >>a = [-10:0.2:10] >>for I = 1:length(a) b(i) = f(a(i)) end ;%画出原函数(a,b) >>c = [-5:1:5] >>for i = 1:length( c) d(i) = f(c(i))
多项式插值法和拉格朗日插值
多项式插值法和拉格朗日插值 教案一多项式插值法和拉格朗日插值 基本内容提要 1 多项式插值法的基本概念 2 插值多项式的存在性与唯一性分析 3 拉格朗日插值多 项式的构造及截断误差 4 截断误差的实用估计式 5 逐次线性插值法教学目的和要求 1 熟练掌握多项式插值法的基本概念 2 理解插值多项式的存在性与唯一性 3 掌握拉 格朗日插值法 4 掌握截断误差的估计方法 5 理解逐次线性插值法的基本思想,掌握Aitken逐次线性插值法 6 掌握运用拉格朗 日插值法处理问题的基本过程教学重点 1 拉格朗日插值基函数及拉格朗日插值多项式的构造 2 拉格朗日插值多项式的截断 误差分析 3 逐次线性插值法的基本思想教学难点 1 插值多项式存在唯一性条件的讨论分析 2 插值误差的分析与估计 3 Aitken逐次线性插值法的计算过程课程类型新知识理论课教学方法 结合提问,以讲授法为主教学过程 问题引入 实际问题中许多变量间的依赖关系往往可用数学中的函数概念刻画,但在多数情况下,这些函数的表达式是未知的,或者函数已知,但形式十分复杂。基于未知函数或复杂函数 的某些已知信息,如何构造这些函数的近似表达式?如何计算这些函数在其它点处的函数值?所构造的近似表达式与真实函数的误差是多少?插值理论与方法就是解决这些问题的 有效工具之一。 §2.1 多项式插值 2.1.1 基本概念 假设f(x)是定义在区间[a,b]上的未知或复杂函数,但已知该函数在点a≤x0 P(xi)=yi,i=0,1,2,L,n,即在给定点xi处,P(x)与f(x)是相吻合的。 (2.1) 把P(x)称为f(x)的插值多项式(函通常把上述x0 数), f(x)称为被插函数。[a,b]称为插值区间,条件(2.1)称为插值条件,并把 求P(x)的过程称为插值法。
Matlab中龙格-库塔(Runge-Kutta)方法原理及实现
Matlab中龙格-库塔(Runge-Kutta)方法原理及实现 龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分,它在积分区间多预计算出几个点的斜率,然后进行加权平均,用做下一点的依据,从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法,如果预先取四个点就是四阶龙格库塔法。一阶常微分方程可以写作:y'=f(x,y),使用差分概念。 (Yn+1-Yn)/h= f(Xn,Yn)推出(近似等于,极限为Yn') Yn+1=Yn+h*f(Xn,Yn) 另外根据微分中值定理,存在0所以,为了更好更准确地把握时间关系,应自己在理解龙格库塔原理的基础上,编写定步长的龙格库塔函数,经过学习其原理,已经完成了一维的龙格库塔函数。 仔细思考之后,发现其实如果是需要解多个微分方程组,可以想象成多个微分方程并行进行求解,时间,步长都是共同的,首先把预定的初始值给每个微分方程的第一步,然后每走一步,对多个微分方程共同求解。想通之后发现,整个过程其实很直观,只是不停的逼近计算罢了。编写的定步长的龙格库塔计算函数: function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称,初始值向量,步长,时间起点,时间终点(参数形式参考了ode45函数) n=floor((b-a)/h);%求步数 x(1)=a;%时间起点 y(:,1)=y0;%赋初值,可以是向量,但是要注意维数 for ii=1:n x(ii+1)=x(ii)+h; k1=ufunc(x(ii),y(:,ii)); k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2); k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2); k4=ufunc(x(ii)+h,y(:,ii)+h*k3); y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6; %按照龙格库塔方法进行数值求解
计算方法上机实验报告——拉格朗日插值问题
计算方法上机实验报告——拉格朗日插值问题 一、方法原理 n次拉格朗日插值多项式为:Ln(x)=y0l0(x)+y1l1(x)+y2l2(x)+…+ynln(x) n=1时,称为线性插值,L1(x)=y0(x-x1)/(x0-x1)+y1(x-x0)/(x1-x0)=y0+(y1-x0)(x-x0)/(x1-x0) n=2时,称为二次插值或抛物线插值,精度相对高些 L2(x)=y0(x-x1)(x-x2)/(x0-x1)/(x0-x2)+y1(x-x0)(x-x2)/(x1-x0)/(x1-x 2)+y2(x-x0)(x-x1)/(x2-x0)/(x2-x1) 二、主要思路 使用线性方程组求系数构造插值公式相对复杂,可改用构造方法来插值。 对节点xi(i=0,1,…,n)中任一点xk(0<=k<=n)作一n次多项式lk(xk),使它在该点上取值为1,而在其余点xi(i=0,1,…,k-1,k+1,…,n)上为0,则插值多项式为Ln(x)=y0l0(x)+y1l1(x)+y2l2(x)+…+ynln(x) 上式表明:n个点xi(i=0,1,…,k-1,k+1,…,n)都是lk(x)的零点。可求得lk 三.计算方法及过程:1.输入节点的个数n 2.输入各个节点的横纵坐标 3.输入插值点 4.调用函数,返回z 函数语句与形参说明 程序源代码如下: 形参与函数类型 参数意义 intn 节点的个数 doublex[n](double*x) 存放n个节点的值 doubley[n](double*y) 存放n个节点相对应的函数值 doublep 指定插值点的值 doublefun() 函数返回一个双精度实型函数值,即插值点p处的近似函数值 #include #include usingnamespacestd; #defineN100 doublefun(double*x,double*y,intn,doublep); voidmain() {inti,n; cout<<"输入节点的个数n:"; cin>>n;
龙格-库塔法MATLAB
1. matlab 新建.m文件,编写龙格-库塔法求解函数 function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称,初始值向量,步长,时间起点,时间终点(参数形式参考了ode45函数) n=floor((b-a)/h); %求步数 x(1)=a;%时间起点 y(:,1)=y0;%赋初值,可以是向量,但是要注意维数 for ii=1:n x(ii+1)=x(ii)+h; k1=ufunc(x(ii),y(:,ii)); k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2); k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2); k4=ufunc(x(ii)+h,y(:,ii)+h*k3); y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6; %按照龙格库塔方法进行数值求解 end 2.另外再新建一个.,m文件,定义要求解的常微分方程函数 function dx=fun1(t,x) dx =zeros(2,1);%初始化列向量 dx(1) =0.08574*x(2)-1.8874*x(1)-20.17; dx(2) =1.8874*x(1)-0.08574*x(2)+115.16; 3,再新建一个.m文件,利用龙格-库塔方法求解常微分方程 [T1,F1]=runge_kutta1(@fun1,[46.30 1296 ],1,0,20); %求解步骤2定义的fun1常微分方程,@fun1是调用其函数指针,从0到20,间隔为1 subplot(122) plot(T1,F1)%自编的龙格库塔函数效果 title('自编的龙格库塔函数') grid 运行步骤3文件即可得到结果,F1为估测值 或者可以调用matlab自带函数ode45求解 方法如下:
数值分析实验一——拉格朗日插值算法报告
拉格朗日插值算法的实现 实验报告 姓名:** 年级:****专业:计算机科学与技术科目:数值分析题目:拉格朗日插值算法的实现 实验时间: 2014年5月27日实验成绩: 实验教师: 一、实验名称:拉格朗日插值算法的实现 二、实验目的: a. 验证拉格朗日插值算法对于不同函数的插值 b. 验证随着插值结点的增多插值曲线的变化情况。 三、实验内容: 拉格朗日插值基函数的一般形式: 也即是: 所以可以得出拉格朗日插值公式的一般形式: 其中, n=1时,称为线性插值,P1(x) = y0*l0(x) + y1*l1(x) n=2时,称为二次插值或抛物插值,精度相对高些,P2(x) = y0*l0(x) + y1*l1(x) + y2*l2(x) 四、程序关键语句描写 double Lagrange(int n,double X[],double Y[],double x) { double result=0; for (int i=0;ifor(int j=0;j #include using namespace std; int main() { double Lagrange(int n,double X[],double Y[],double x); //插值函数double x;//要求插值的x的值 double result;//插值的结果 char a='n'; double X[20],Y[20]; do { cout<<"请输入插值次数n的值:"<>n; cout<<"请输入插值点对应的值及函数值(xi,yi):"<>X[k]>>Y[k]; } cout<<"请输入要求值x的值:"<>x; result=Lagrange(n,X,Y,x); cout<<"由拉格朗日插值法得出结果:"<>a; }while(a=='yes'); return 0; }
拉格朗日插值龙格现象的matlab实现
拉格朗日插值龙格现象的MATLAB 实现 姓名:袁宽 学号:2 专业:电气工程 题目:对于函数2 11)(x x f +=,55≤≤-x 进行拉格朗日插值。10=n ,按等距节点求分段线性插值,把)(x f 和插值多项式的曲线画在同一张图上进行比较。观察Lagrange 插值及数值积分中的分段性插值。 f.m : function f= f( x ) f=1./(1+x.^2); end Lagrange.m function y=Lagrange(x0,y0,x); n=length(x0); m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=k p=p*(z-x0(j))/(x0(k)-x0(j)); end end
s=p*y0(k)+s; end y(i)=s; End 拉格朗日插值的曲线: x=[-5:1:5]; y=1./(1+x.^2); x0=[-5:0.001:5]; y0=Lagrange(x,y,x0); y1=1./(1+x0.^2); plot(x0,y0,'b') hold on plot(x0,y1,'r') 使用龙格现象观察分段性插值的曲线:syms y x lx; y=1/(1+x^2); x0=-5:1:5; y0=zeros(1,length(x0)); for i=1:11 x=x0(i);
y0(i)=eval(y); end M = -5:0.01:5; y1 = zeros(1,length(M)); n = 1; for i=2:11 for x=-5:0.01:5 if x=x0(i-1) lx(1)=y0(i-1)*(x-x0(i))/(x0(i-1)-x0(i)); lx(2)=y0(i)*(x-x0(i-1))/(x0(i)-x0(i-1)); y1(n) = lx(1)+lx(2); n = n+1; end end end ezplot(y,[-5,5]) hold on A =-5:0.01:5; plot(A,y1,'r'); 分别运行两个文件后的得到两个图形曲线:
龙格现象实验报告1
数值分析实验报告 实验名称:观察龙格(Runge)现象实验 班级:12级信息与计算科学(1)班 姓名: 学号:33 16 17 59 实验日期: 2014.10.11 周次: 6 实验地点: A14-504
多种插值对比 function y=fun(x); y=5./(1+x.^2); end hours=-5:5; h=-5:0.1:5; temps=fun(hours); t1=interp1(hours,temps,h,'spline'); %(?±?óê?3?êy?Y??ê?oü?àμ?) t2=interp1(hours,temps,h, 'nearest'); t3=interp1(hours,temps,h, 'linear'); t4=interp1(hours,temps,h, 'cubic'); subplot(2,2,1); plot(hours,temps, ' bo',h,t1, 'r'); title('spline '); %×÷í? subplot(2,2,2); plot(hours,temps, 'bo',h,t2, 'r'); title(' nearest'); subplot(2,2,3); plot(hours,temps, 'bo',h,t3, 'r'); title('linear'); subplot(2,2,4); plot(hours,temps, 'bo',h,t4, 'r', h,t1, 'g'); title('cubic-spline'); -5 05 spline -5 05 nearest -5 05 linear -5 05 cubic-spline
龙格库塔法求微分方程2
《MATLAB 程序设计实践》课程考核 一、编程实现“四阶龙格-库塔(R-K )方法求常微分方程”,并举一 例应用之。 【实例】采用龙格-库塔法求微分方程: ?? ?==+-=0 , 0)(1 '00 x x y y y 1、算法说明: 在龙格-库塔法中,四阶龙格-库塔法的局部截断误差约为o(h5),被广泛应用于解微分方程的初值问题。其算法公式为: )22(6 3211k k k h y y n n +++=+ 其中: ?????????++=++=++ ==) ,() 21 ,21()21 ,21() ,(34 23121hk y h x f k hk y h x f k hk y h x f k y x f k n n n n n n n n 2、流程图: 2.1、四阶龙格-库塔(R-K )方法流程图:
2.2、实例求解流程图:
3、源程序代码 3.1、四阶龙格-库塔(R-K)方法源程序: function [x,y] = MyRunge_Kutta(fun,x0,xt,y0,PointNum,varargin) %Runge-Kutta 方法解微分方程形为 y'(t)=f(x,y(x)) %此程序可解高阶的微分方程。只要将其形式写为上述微分方程的向量形式 %函数 f(x,y): fun %自变量的初值和终值:x0, xt %y0表示函数在x0处的值,输入初值为列向量形式 %自变量在[x0,xt]上取的点数:PointNum %varargin为可输入项,可传适当参数给函数f(x,y) %x:所取的点的x值 %y:对应点上的函数值 if nargin<4 | PointNum<=0 PointNum=100; end if nargin<3 y0=0; end y(1,:)=y0(:)'; %初值存为行向量形式h=(xt-x0)/(PointNum-1); %计算步长 x=x0+[0:(PointNum-1)]'*h; %得x向量值 for k=1:(PointNum)%迭代计算 f1=h*feval(fun,x(k),y(k,:),varargin{:}); f1=f1(:)'; %得公式k1 f2=h*feval(fun,x(k)+h/2,y(k,:)+f1/2,varargin{:}); f2=f2(:)'; %得公式k2 f3=h*feval(fun,x(k)+h/2,y(k,:)+f2/2,varargin{:}); f3=f3(:)'; %得公式k3 f4=h*feval(fun,x(k)+h,y(k,:)+f3,varargin{:}); f4=f4(:)'; %得公式k4 y(k+1,:)=y(k,:)+(f1+2*(f2+f3)+f4)/6; %得y(n+1) end 3.2、实例求解源程序: %运行四阶R-K法