数值分析龙格现象matlab代码分享
数值分析matlab作业龙格库塔欧拉方法解二阶微分方程
Matlab 应用使用Euler 和Rungkutta 方法解臂状摆的能量方程背景 单摆是常见的物理模型，为了得到摆角θ的关于时间的函数，来描述单摆运动。
由角动量定理我们知道εJ M =化简得到0sin 22=+θθlgdt d在一般的应用和计算中，只考虑摆角在5度以内的小摆动，因为可以吧简化为θ，这样比较容易解。
实际上这是一个解二阶常微分方程的问题。
在这里的单摆是一种特别的单摆，具有均匀的质量M 分布在长为2的臂状摆上， 使用能量法建立方程WT =hmg ∆=2J 21ω 化简得到θθcos 35499.722=dtd重力加速度取9.806651使用欧拉法令dxdy z =，这样降阶就把二阶常微分方程转化为一阶微分方程组，再利用向前Euler 方法数值求解。
y(i+1)=y(i)+h*z(i);z(i+1)=z(i)+h*7.35499*cos(y(i)); y(0)=0 z(0)=0精度随着h 的减小而更高，因为向前欧拉方法的整体截断误差与h 同阶，（因为是用了泰勒公式）所以欧拉方法的稳定区域并不大。
2.RK4-四阶龙格库塔方法使用四级四阶经典显式Rungkutta 公式稳定性很好，RK4法是四阶方法，每步的误差是h5阶，而总积累误差为h4阶。
所以比欧拉稳定。
运行第三个程序：在一幅图中显示欧拉法和RK4法，随着截断误差的积累，欧拉法产生了较大的误差h=0.01h=0.0001，仍然是开始较为稳定，逐渐误差变大总结：RK4是很好的方法，很稳定，而且四阶是很常用的方法，因为到五阶的时候精度并没有相应提升。
通过这两种方法计算出角度峰值y=3.141593，周期是1.777510。
三个程序欧拉法clear;clch=0.00001;a=0;b=25;x=a:h:b;y(1)=0;z(1)=0;for i=1:length(x)-1 % 欧拉y(i+1)=y(i)+h*z(i);z(i+1)=z(i)+h*7.35499*cos(y(i));endplot(x,y,'r*');xlabel('时间');ylabel('角度');A=[x,y];%y(find(y==max(y)))%Num=(find(y==max(y)))[y,T]=max(y);fprintf('角度峰值等于%d',y) %角度的峰值也就是πfprintf('\n')fprintf('周期等于%d',T*h) %周期legend('欧拉');龙格库塔方法先定义函数rightf_sys1.mfunction w=rightf_sys1(x,y,z)w=7.35499*cos(y);clear;clc;%set(0,'RecursionLimit',500)h=0.01;a=0;b=25;x=a:h:b;RK_y(1)=0; %初值RK_z(1)=0; %初值for i=1:length(x)-1K1=RK_z(i);L1=rightf_sys1(x(i),RK_y(i),RK_z(i)); % K1 and L1 K2=RK_z(i)+0.5*h*L1;L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1);K3=RK_z(i)+0.5*h*L2;L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2);K4=RK_z(i)+h*L3;L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); % K4 and L4RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);endplot(x,RK_y,'b+');xlabel('Variable x');ylabel('Variable y');A=[x,RK_y];[y,T]=max(RK_y);legend('RK4方法');fprintf('角度峰值等于%d',y) %角度的峰值也就是πfprintf('\n')fprintf('周期等于%d',T*h) %周期两个方法在一起对比使用跟上一个相同的函数rightf_sys1.mclear;clc; %清屏h=0.0001;a=0;b=25;x=a:h:b;Euler_y(1)=0;Euler_z(1)=0; %欧拉的初值RK_y(1)=0;RK_z(1)=0; %龙格库塔初值for i=1:length(x)-1%先是欧拉法Euler_y(i+1)=Euler_y(i)+h*Euler_z(i);Euler_z(i+1)=Euler_z(i)+h*7.35499*cos(Euler_y(i));%龙格库塔K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i)); % K1 and L1 K2=RK_z(i)+0.5*h*L1;L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1);% K2 and L2K3=RK_z(i)+0.5*h*L2;L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2);% K3 and L3K4=RK_z(i)+h*L3;L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); % K4 and L4RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);endplot(x,Euler_y,'r-',x,RK_y,'b-');[y,T]=max(RK_y);fprintf('角度峰值等于%d',y) %角度的峰值也就是πfprintf('\n')fprintf('周期等于%d',T*h) %周期xlabel('时间');ylabel('角度');legend('欧拉','RK4');文- 汉语汉字编辑词条文，wen，从玄从爻。
数值分析MATLAB编程-推荐下载
题目三：设������(������) = ������������，在[-1,1]上用 Legendre 多项式作������(������)的 3 次最佳平方逼近多项式。 >> a=-1;b=1; >> n=3; >> syms x; >> fx=exp(x); >> exp(x); >> for k=3 F=x^k*fx; d(k+1)=int(F,a,b); end >> d=reshape(d,n+1,1); >> H=hilb(n+1); >Fra bibliotek a=H\d a=
对全部高中资料试卷电气设备，在安装过程中以及安装结束后进行高中资料试卷调整试验；通电检查所有设备高中资料电试力卷保相护互装作置用调与试相技互术通关，1系电过，力管根保线据护敷生高设产中技工资术艺料0不高试仅中卷可资配以料置解试技决卷术吊要是顶求指层，机配对组置电在不气进规设行范备继高进电中行保资空护料载高试与中卷带资问负料题荷试2下卷2，高总而中体且资配可料置保试时障卷，各调需类控要管试在路验最习；大题对限到设度位备内。进来在行确管调保路整机敷使组设其高过在中程正资1常料中工试，况卷要下安加与全强过，看度并22工且22作尽22下可22都能22可地护以缩1关正小于常故管工障路作高高；中中对资资于料料继试试电卷卷保破连护坏接进范管行围口整，处核或理对者高定对中值某资，些料审异试核常卷与高弯校中扁对资度图料固纸试定，卷盒编工位写况置复进.杂行保设自护备动层与处防装理腐置，跨高尤接中其地资要线料避弯试免曲卷错半调误径试高标方中高案资等，料，编试要5写、卷求重电保技要气护术设设装交备备置底4高调、动。中试电作管资高气，线料中课并敷3试资件且、设卷料中拒管技试试调绝路术验卷试动敷中方技作设包案术，技含以来术线及避槽系免、统不管启必架动要等方高多案中项；资方对料式整试，套卷为启突解动然决过停高程机中中。语高因文中此电资，气料电课试力件卷高中电中管气资壁设料薄备试、进卷接行保口调护不试装严工置等作调问并试题且技，进术合行，理过要利关求用运电管行力线高保敷中护设资装技料置术试做。卷到线技准缆术确敷指灵设导活原。。则对对：于于在调差分试动线过保盒程护处中装，高置当中高不资中同料资电试料压卷试回技卷路术调交问试叉题技时，术，作是应为指采调发用试电金人机属员一隔，变板需压进要器行在组隔事在开前发处掌生理握内；图部同纸故一资障线料时槽、，内设需，备要强制进电造行回厂外路家部须出电同具源时高高切中中断资资习料料题试试电卷卷源试切，验除线报从缆告而敷与采设相用完关高毕技中，术资要资料进料试行，卷检并主查且要和了保检解护测现装处场置理设。备高中资料试卷布置情况与有关高中资料试卷电气系统接线等情况，然后根据规范与规程规定，制定设备调试高中资料试卷方案。
数值分析matlab程序
数值分析（matlab程序）曹德欣曹璎珞第一章绪论数值稳定性程序，计算P11 试验题一积分function try_stableglobal nglobal aa=input('a=');N = 20;I0 = log(1+a)-log(a);I = zeros(N,1);I(1) = -a*I0+1;for k = 2:NI(k) = - a*I(k-1)+1/k;endII = zeros(N,1);if a>=N/(N+1)II(N) = (1+2*a)/(2*a*(a+1)*(N+1));elseII(N) =(1/((a+1)*(N+1))+1/N)/2;endfor k = N:-1:2II(k-1) = ( - II(k) +1/k) / a;endIII = zeros(N,1);for k = 1:Nn = k;III(k) = quadl(@f,0,1);endclcfprintf('\n 算法1结果算法2结果精确值') for k = 1:N,fprintf('\nI(%2.0f) %17.7f %17.7f %17.7f',k,I(k),II(k),III(k)) endfunction y = f(x)global nglobal ay = x.^n./(a+x);return第二章非线性方程求解下面均以方程y=x^4+2*x^2-x-3为例：1、二分法function y=erfen(a,b,esp)format longif nargin<3 esp=1.0e-4;endif fun(a)*fun(b)<0n=1;c=(a+b)/2;while c>espif fun(a)*fun(c)<0b=c;c=(a+b)/2;elseif fun(c)*fun(b)<0a=c;c=(a+b)/2;else y=c; esp=10000;endn=n+1;endy=c;elseif fun(a)==0y=a;elseif fun(b)==0y=b;else disp('these,nay not be a root in the intercal')endnfunction y=fun(x)y=x^4+2*x^2-x-3;2、牛顿法function y=newton(x0)x1=x0-fun(x0)/dfun(x0);n=1;while (abs(x1-x0)>=1.0e-4) & (n<=100000000)x0=x1;x1=x0-fun(x0)/dfun(x0);n=n+1;endy=x1nfunction y=fun(x)y=x^4+2*x^2-x-3; 3、割线法function y=gexian(x0,x1)x2=x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0)); %根据初始XO 和X1求X2 n=1;while (abs(x1-x0)>=1.0e-4) & (n<=100000000) %判断两个条件截止 x0=x1; %将x1赋给x0 x1=x2; %将x2赋给x1 x2=x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0)); %迭代运算 n=n+1; end y=x2 nfunction y=fun(x)y=x^4+2*x^2-x-3;第四章题目：推导外推样条公式：⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--------1232123211223322~~~~~22~n n n n n n n n d d d d M M M Mδμλμλμλδ，并编写程序与Matlab 的Spline 函数结果进行对比，最后调用追赶法解方程组。
数值分析matlab程序
数值分析（matlab程序） 曹德欣 曹璎珞 第一章 绪论 数值稳定性程序，计算P11 试验题一积分 function try_stable global n global a a=input('a='); N = 20; I0 = log(1+a)-log(a); I = zeros(N,1); I(1) = -a*I0+1; for k = 2:N I(k) = - a*I(k-1)+1/k; end II = zeros(N,1); if a>=N/(N+1) II(N) = (1+2*a)/(2*a*(a+1)*(N+1)); else II(N) =(1/((a+1)*(N+1))+1/N)/2; end for k = N:-1:2 II(k-1) = ( - II(k) +1/k) / a; end III = zeros(N,1); for k = 1:N n = k; III(k) = quadl(@f,0,1); end clc fprintf('\n 算法1结果 算法2结果 精确值') for k = 1:N, fprintf('\nI(%2.0f) %17.7f %17.7f %17.7f',k,I(k),II(k),III(k)) end function y = f(x) global n global a y = x.^n./(a+x); return 第二章 非线性方程求解 下面均以方程y=x^4+2*x^2-x-3为例： 1、二分法 function y=erfen(a,b,esp) format long if nargin<3 esp=1.0e-4; end if fun(a)*fun(b)<0 n=1; c=(a+b)/2; while c>esp if fun(a)*fun(c)<0 b=c; c=(a+b)/2; elseif fun(c)*fun(b)<0 a=c; c=(a+b)/2; else y=c; esp=10000; end n=n+1; end y=c; elseif fun(a)==0 y=a; elseif fun(b)==0 y=b; else disp('these,nay not be a root in the intercal') end n function y=fun(x) y=x^4+2*x^2-x-3; 2、牛顿法 function y=newton(x0) x1=x0-fun(x0)/dfun(x0); n=1; while (abs(x1-x0)>=1.0e-4) & (n<=100000000) x0=x1; x1=x0-fun(x0)/dfun(x0); n=n+1; end y=x1 n function y=fun(x) y=x^4+2*x^2-x-3; 3、割线法 function y=gexian(x0,x1) x2=x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0)); %根据初始XO和X1求X2 n=1; while (abs(x1-x0)>=1.0e-4) & (n<=100000000) %判断两个条件截止 x0=x1; %将x1赋给x0 x1=x2; %将x2赋给x1 x2=x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0)); %迭代运算 n=n+1; end y=x2 n function y=fun(x) y=x^4+2*x^2-x-3;
MATLAB常微分方程数值解——欧拉法、改进的欧拉法与四阶龙格库塔方法
MATLAB常微分⽅程数值解——欧拉法、改进的欧拉法与四阶龙格库塔⽅法MATLAB常微分⽅程数值解作者：凯鲁嘎吉 - 博客园1.⼀阶常微分⽅程初值问题2.欧拉法3.改进的欧拉法4.四阶龙格库塔⽅法5.例题⽤欧拉法，改进的欧拉法及4阶经典Runge-Kutta⽅法在不同步长下计算初值问题。
步长分别为0.2,0.4,1.0.matlab程序：function z=f(x,y)z=-y*(1+x*y);function R_K(h)%欧拉法y=1;fprintf('欧拉法：x=%f, y=%f\n',0,1);for i=1:1/hx=(i-1)*h;K=f(x,y);y=y+h*K;fprintf('欧拉法：x=%f, y=%f\n',x+h,y);endfprintf('\n');%改进的欧拉法y=1;fprintf('改进的欧拉法：x=%f, y=%f\n',0,1);for i=1:1/hx=(i-1)*h;K1=f(x,y);K2=f(x+h,y+h*K1);y=y+(h/2)*(K1+K2);fprintf('改进的欧拉法：x=%f, y=%f\n',x+h,y);endfprintf('\n');%龙格库塔⽅法y=1;fprintf('龙格库塔法：x=%f, y=%f\n',0,1);for i=1:1/hx=(i-1)*h;K1=f(x,y);K2=f(x+h/2,y+(h/2)*K1);K3=f(x+h/2,y+(h/2)*K2);K4=f(x+h,y+h*K3);y=y+(h/6)*(K1+2*K2+2*K3+K4);fprintf('龙格库塔法：x=%f, y=%f\n',x+h,y);end结果：>> R_K(0.2)欧拉法：x=0.000000, y=1.000000欧拉法：x=0.200000, y=0.800000欧拉法：x=0.400000, y=0.614400欧拉法：x=0.600000, y=0.461321欧拉法：x=0.800000, y=0.343519欧拉法：x=1.000000, y=0.255934改进的欧拉法：x=0.000000, y=1.000000改进的欧拉法：x=0.200000, y=0.807200改进的欧拉法：x=0.400000, y=0.636118改进的欧拉法：x=0.600000, y=0.495044改进的欧拉法：x=0.800000, y=0.383419改进的欧拉法：x=1.000000, y=0.296974龙格库塔法：x=0.000000, y=1.000000龙格库塔法：x=0.200000, y=0.804636龙格库塔法：x=0.400000, y=0.631465龙格库塔法：x=0.600000, y=0.489198龙格库塔法：x=0.800000, y=0.377225龙格库塔法：x=1.000000, y=0.291009>> R_K(0.4)欧拉法：x=0.000000, y=1.000000欧拉法：x=0.400000, y=0.600000欧拉法：x=0.800000, y=0.302400改进的欧拉法：x=0.000000, y=1.000000改进的欧拉法：x=0.400000, y=0.651200改进的欧拉法：x=0.800000, y=0.405782龙格库塔法：x=0.000000, y=1.000000龙格库塔法：x=0.400000, y=0.631625龙格库塔法：x=0.800000, y=0.377556>> R_K(1)欧拉法：x=0.000000, y=1.000000欧拉法：x=1.000000, y=0.000000改进的欧拉法：x=0.000000, y=1.000000改进的欧拉法：x=1.000000, y=0.500000龙格库塔法：x=0.000000, y=1.000000龙格库塔法：x=1.000000, y=0.303395注意：在步长h为0.4时，要将for i=1:1/h改为for i=1:0.8/h。
拉格朗日插值、牛顿插值的matlab代码
实验五多项式插值逼近信息与计算科学金融崔振威201002034031一、实验目的：拉格朗日插值和牛顿插值的数值实现二、实验内容：p171.1、p178.1、龙格现象数值实现三、实验要求：1、根据所给题目构造相应的插值多项式，2、编程实现两类插值多项式的计算3、试分析多项式插值造成龙格现象的原因主程序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+1v=1;for j=1:n+1if k~=jv=conv(v,poly(x(j)))/(x(k)-x(j)) %对多项式做卷积运算endendl(k,:)=v;endc=y*l;牛顿插值多项式主程序function [p2,z]=newTon(x,y,t)%输入参数中x,y为元素个数相等的向量%t为插入的定点%p2为所求得的牛顿插值多项式%z为利用多项式所得的t的函数值。
n=length(x);chaS(1)=y(1);for i=2:nx1=x;y1=y;x1(i+1:n)=[];y1(i+1:n)=[];n1=length(x1);s1=0;for j=1:n1t1=1;for k=1:n1if k==j %如果相等则跳出循环continue;elset1=t1*(x1(j)-x1(k));endends1=s1+y1(j)/t1;endchaS(i)=s1;endb(1,:)=[zeros(1,n-1) chaS(1)];cl=cell(1,n-1); %cell定义了一个矩阵for i=2:nu1=1;for j=1:i-1u1=conv(u1,[1 -x(j)]); %conv()用于多项式乘法、矩阵乘法cl{i-1}=u1;endcl{i-1}=chaS(i)*cl{i-1};b(i,:)=[zeros(1,n-i),cl{i-1}];endp2=b(1,:);for j=2:np2=p2+b(j,:);endif length(t)==1rm=0;for i=1:nrm=rm+p2(i)*t^(n-i);endz=rm;elsek1=length(t);rm=zeros(1,k1);for j=1:k1for i=1:nrm(j)=rm(j)+p2(i)*t(j)^(n-i);endz=rm;endendplot(t,z,'y',x,y,'*r') %输出牛顿插值多项式的函数图p171.1(a)、f(x)=e x解：在matlab窗口中输入：>> x=[0 0.2 0.4 0.6 0.8 1];>> y=[exp(0) exp(0.2) exp(0.4) exp(0.6) exp(0.8) exp(1)]y =1.0000 1.2214 1.4918 1.82212.2255 2.7183>> [c,l]=lagran(x,y)可以得出输出结果为：c =0.0139 0.0349 0.1704 0.4991 1.0001 1.0000l =-26.0417 78.1250 -88.5417 46.8750 -11.4167 1.0000130.2083 -364.5833 369.7917 -160.4167 25.0000 0-260.4167 677.0833 -614.5833 222.9167 -25.0000 0260.4167 -625.0000 510.4167 -162.5000 16.6667 0-130.2083 286.4583 -213.5417 63.5417 -6.2500 026.0417 -52.0833 36.4583 -10.4167 1.0000 0由输出结果可以的出：P(x)的系数分别为：a0=0.0139 a1=0.0349 a2=0.1704 a3=0.4991 a4=1.0001 a5=1.0000(b)、f(x)=sin(x)解：在matlab窗口中输入：>> x=[0 0.2 0.4 0.6 0.8 1];>> y=[sin(0) sin(0.2) sin(0.4) sin(0.6) sin(0.8) sin(1)];>> [c,l]=lagran(x,y)可以得出输出结果为：c =0.0073 0.0016 -0.1676 0.0002 1.0000 0l =-26.0417 78.1250 -88.5417 46.8750 -11.4167 1.0000130.2083 -364.5833 369.7917 -160.4167 25.0000 0-260.4167 677.0833 -614.5833 222.9167 -25.0000 0260.4167 -625.0000 510.4167 -162.5000 16.6667 0-130.2083 286.4583 -213.5417 63.5417 -6.2500 026.0417 -52.0833 36.4583 -10.4167 1.0000 0由输出结果可以的出：P(x)的系数分别为：a0=0.0073 a1=0.0016 a2=-0.1676 a3=0.0002 a4=1.0000 a5=0(c)、f(x)=(x+1)x+1解：在matlab窗口中输入：>> x=[0 0.2 0.4 0.6 0.8 1];>> y=[1 1.2^1.2 1.4^1.4 1.6^1.6 1.8^1.8 2^2];>> [c,l]=lagran(x,y)可以得出输出结果为：c =0.3945 -0.0717 0.7304 0.9415 1.0052 1.0000l =-26.0417 78.1250 -88.5417 46.8750 -11.4167 1.0000130.2083 -364.5833 369.7917 -160.4167 25.0000 0-260.4167 677.0833 -614.5833 222.9167 -25.0000 0260.4167 -625.0000 510.4167 -162.5000 16.6667 0-130.2083 286.4583 -213.5417 63.5417 -6.2500 026.0417 -52.0833 36.4583 -10.4167 1.0000 0由输出结果可以的出：P(x)的系数分别为：a0=0.3945 a1=-0.0717 a2=0.7304 a3=0.9415 a4=1.0052 a5=1.0000P178.12、a0=5 a1=-2 a2=0.5 a3=-0.1 a4=0.003x0=0 x1=1 x2=2 x3=3 c=2.5解：在matlab窗口中输入：>> x=[5 -2 0.5 -0.1];>> y=[0 1 2 3];>> t=0:0.1:2.5;>> [u,v]=newTon(x,y,t)可得出输出结果：u =0.1896 -0.7843 -1.3928 2.8688v =2.8688 2.7218 2.5603 2.3855 2.1983 2.0000 1.7917 1.5745 1.3497 1.1182 0.8813 0.6401 0.3957 0.1493 -0.0980 -0.3451 -0.5908 -0.8340 -1.0735 -1.3082 -1.5370 -1.7588 -1.9723 -2.1765 -2.3702 -2.5523由此可以求出牛顿多项式为：f(x)=0.1896x^3--0.7843^x2--1.3928x+2.8688输出的图为：结果分析:利用牛顿插值多项式的函数，通过调用函数可以求得牛顿多项式与给定的点的值，并通过matlab做出函数图像。
数值积分龙贝格matlab
《数值分析》课程实验报告一、实验目的1、进一步熟悉向量矩阵的运算；2、掌握龙贝格（Romberg ）算法，并能用高级程序语言MATLAB 编写实现此算法的程序；3、进而加深对龙贝格（Romberg ）算法的理解。
二、实验内容1. 使用Romberg 积分，对于计算下列⎰+4802)cos (1dx x 各近似值a.确定1,51,41,31,21,1,,,,R R R R Rb.确定5,54,43,32,2,,,R R R Rc.6,65,64,63,62,61,6,,,,,R R R R R Rd.确定10,109,98,87,7,,,R R R R三、实验步骤1. 编写程序龙贝格积分方法如下：n=5;a=0;b=48;h(1,1)=b-a;fa=sqrt(1+(cos(a))^2);fb=sqrt(1+(cos(b))^2);r(1,1)=h(1,1)/2*(fa+fb);disp('R11，R21，R31，R41，R51分别为');disp(r(1,1));for i=2:nh(i,1)=(b-a)/(2^(i-1));sum=0;for k=1:2^(i-2)x=a+(2*k-1)*h(i,1);sum=sum+sqrt(1+(cos(x)).^2);endr(i,1)=0.5*(r(i-1,1)+h(i-1,1)*sum);disp(r(i,1));enddisp('R22，R33，R44，R55分别为');for k=2:nfor j=2:kr(k,j)=r(k,j-1)+(r(k,j-1)-r(k-1,j-1))/(4^(j-1)-1);enddisp(r(k,k));enddisp('R61，R62，R63，R64，R65，R66分别为');n=6;for i=2:nh(i,1)=(b-a)/(2^(i-1));sum=0;for k=1:2^(i-2)x=a+(2*k-1)*h(i,1);sum=sum+sqrt(1+(cos(x)).^2);endr(i,1)=0.5*(r(i-1,1)+h(i-1,1)*sum);endfor k=2:nfor j=2:kr(k,j)=r(k,j-1)+(r(k,j-1)-r(k-1,j-1))/(4^(j-1)-1);endendfor i=1:ndisp(r(6,i));enddisp('R77，R88，R99，R10,10分别为');n=10;for i=2:nh(i,1)=(b-a)/(2^(i-1));sum=0;for k=1:2^(i-2)x=a+(2*k-1)*h(i,1);sum=sum+sqrt(1+(cos(x)).^2);endr(i,1)=0.5*(r(i-1,1)+h(i-1,1)*sum);endfor k=2:nfor j=2:kr(k,j)=r(k,j-1)+(r(k,j-1)-r(k-1,j-1))/(4^(j-1)-1);endendfor i=7:10disp(r(i,i));end运行结果如下：R11，R21，R31，R41，R51分别为62.437457.288656.443856.263156.2188R22，R33，R44，R55分别为55.572356.201556.205656.2041R61，R62，R63，R64，R65，R66分别为58.362759.077359.268959.317559.329759.3328R77，R88，R99，R10,10分别为58.422158.470758.470558.4705四、实验小结在这次编程中我学到了很多东西，把程序写进软件中也出现了很多错误，细节问题使我们必须注意的，自己有了很多的收获，自己进一步理解和学习了Matlab软件。
matlab迭龙格库塔法解常微分方程
一、介绍迭龙格-库塔法（Runge-Kutta method）是一种数值求解常微分方程（ODE）的常用方法。
它是由卡尔·迭龙格（Carl Runge）和马丁·威尔黑尔姆·库塔（Wilhelm Kutta）在20世纪初提出的，该方法以两位数值分析家的名字来命名。
二、简单描述迭龙格-库塔法是通过数值逼近的方式，来计算常微分方程的近似解。
它是一种显式求解方法，适用于解非线性常微分方程和具有较大阶数的常微分方程。
三、数学原理迭龙格-库塔法主要是通过将微分方程转化为差分方程，利用数值解的方式来逼近微分方程的解。
它是一种显式方法，通过不断迭代得到下一个时间步的近似解。
四、matlab中的应用在matlab中，可以使用ode45函数来调用迭龙格-库塔法求解常微分方程。
ode45函数是matlab中集成的一个函数，通过调用ode45函数，可以直接求解常微分方程的数值解。
五、实例演示下面通过一个简单的例子来演示如何使用matlab中的ode45函数来求解常微分方程。
我们考虑一个简单的一阶常微分方程：dy/dt = -y初始条件为y(0) = 1。
在matlab中，可以通过以下代码来求解该微分方程：```定义微分方程的函数function dydt = myode(t, y)dydt = -y;调用ode45函数求解[t, y] = ode45(myode, [0, 5], 1);plot(t, y);```运行以上代码，即可得到微分方程的数值解，并通过绘图来展示解的变化。
六、总结迭龙格-库塔法是一种常用的数值解常微分方程的方法，它在matlab中有较为方便的调用方式。
通过ode45函数，可以快速求解常微分方程的数值解，并通过绘图来展示结果。
希望本篇文章对读者有所帮助，谢谢阅读。
七、应用场景和优势在实际应用中，迭龙格-库塔法广泛应用于各种科学和工程领域，如物理学、化学、生物学、经济学等。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
1
龙格现象(
Runge phenomenon)
高次插值的病态性质
1. 先建立一个n+1个插值节点的拉格朗日插值多项式
function langrange= langrange( x,n )
langrange=0;
xx=linspace(-5,5,n+1);
for i=1:n+1
lix=1;
for j=1:n+1
if j~=i
lix=lix.*((x-xx(j))./(xx(i)-xx(j)));
end
end
langrange=fun(xx(i)).*lix+langrange;
end
end
2. 再建立一个龙格函数
function f= fun( x )
f=1./(1+x.^2);
3.在同一坐标系中画出龙格函数和拉格朗日插值多项式的图像
function runge_phen(n)
% n为Lagrange插值节点的个数
x=linspace(-5,5,100);
plot(x,fun(x),'r+',x,langrange(x,n),'b*');
2
4.在Matlab命令窗口运行如下命令
runge_phen(10)