matlab龙贝格算法
1. MATLAB用龙贝格积分算法计算f=x^2在0到1上的积分
f = @(x) x.*x % f是被积函数
a = 0;b = 1; % a,b分别是积分的上下限
n = 5; % n+1是T数表的列数
delta = 10^(-8); % 允许误差
M=1;
h=b-a;
err=1;
J=0;
R=zeros(4,4);
R(1,1)=h*(feval(f,a)+feval(f,b))/2;
while ((err>delta)&&(J<n))||(J<4)
J=J+1;
h=h/2;
s=0;
for p=1:M
x=a+h*(2*p-1);
s=s+feval(f,x);
end
R(J+1,1)=R(J,1)/2+h*s;
M=2*M;
for K=1:J
R(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K))/(4^K-1);
end
err=abs(R(J,J)-R(J+1,K+1));
end
quad11=R(J+1,J+1);
quad11 % 所求积分值
R % R是T数表
2. matlab编程 使用龙贝格积分公式 计算定积分 从0-1上的 ln(x)/x 的积分。
因为函数 ln(1+x)/x在[0,1]上是反常积分。
3. 用matlab龙贝格算法
主程序
clc;
clear;
format long e;%精确显示
y=18;
T=zeros(y,y);%定义长度为y的矩阵T
for i=0:y-1
T(i+1)=i;%对矩阵T的第一行赋初值
end
a=0;
b=1;
T(1,2)=(b-a)*(1+f(b))/2;%算出T1
for l=2:y
sum=0;
for j=1:2^(l-1)
sum=sum+f(a+(2*j-1)*(b-a)/2^l);
end
T(l,2)=T(l-1,2)/2+((b-a)/2^l)*sum; %算出梯形序列
if (T(l,2)-T(l-1,2))<0.0000005; %判断是否符合精度
break;
end
end
h=1;
for m=1:y-2
for k=h:(y+h-m-1)
T(k+1,m+2)=((4^m)*T(k+1,m+1)-T(k,m+1))/(4^m-1);
%按公式算出接下去的序列
if (T(k+1,m+2)-T(k,m+2))<0.0000005;
break;
end
end
h=h+1;
end
4. Matlab/simulink中,什么叫Oder45和Ode23bt算法
ode45是基于四点法和五点法的解微分方程数值解的方法,ode23等也一样,都是基于已知点“预测”下一个点的函数值的方法,不同的算法“预测”的方法不一样。比较着名的“预测”方法有欧拉法,改进的欧拉法,龙格库塔法,多点法等。在matlab一般使用中这些方法的差别不大,可以不予理会,会用一个即可,推荐ODE4。
5. 用matlab 龙贝格求积分
我试了下 把你的程序保存为Romberg.m
在工作区输入f=@(x) 1/(x+1)
a=0
b=1
eps=10^(-4)
Romberg(f,a,b,eps)
错误提示是你程序第13行的变量 s没有定义 应该是大写吧
6. 用matlab程序实现:用romberg方法计算积分sinx/x上限为1下限为0的近似值,要求误差不超过0.5*10……-6
a=0;
b=1;
n=1;
k=1;
e=0.5*10^(-6);
T=zeros(10,10);
h=(b-a)/2;
t=h*(1+sin(1));
T(1,1)=t;
for j=1:9
F=0;
for i=1:n
F=F+sin(a+(2*i-1)*h)/(a+(2*i-1)*h);
end
for i=1:k
if T(1,i+1)==0
T(1,i+1)=T(1,i)/2+h*F;
end
end
for m=1:k
if T(m+1,k-m+1)==0
T(m+1,k-m+1)=(4^m*T(m,k-m+2)-T(m,k-m+1))/(4^m-1);
end
end
if abs(T(m+1,1)-T(m,1))<e
I=T(m+1,1);
break
else
h=h/2;
n=2*n;
k=k+1;
end
end
最后的 I 就是函数值,K就是满足精度时的M值
7. 将区间[0 1]等分,并将每一个等分区间进行积分的matlab程序
下面是使用龙贝格算法求积分的matlab程序代码
clear
clc
format long
f='4/(1+x^2)'; %这是被积函数
x='x'; %这是被积自变量
a=0; %这是积分下限
b=1; %这是积分上限
e=1e-5; %这是积分误差限制
%以下是龙贝格积分算法,是目前最为成熟的积分算法,具有收敛速度快,精度可以自定义的优点
% I为积分的估计值
% n为迭代次数,2^(n-1)是等分区间的份数
T(1,1)=(b-a)/2*(subs(f,x,a)+subs(f,x,b));
T=double(T);
n=2;
h=b-a;
T(2,1)=T(1,1)/2+h/2*double(subs(f,x,a+h/2));
T(2,2)=4/3*T(2,1)-1/3*T(1,1);
d=T(2,2)-T(1,1);
while d>e
n=n+1;
h=h/2;
T(n,1)=T(n-1,1)/2;
for i=1:2^(n-2)
T(n,1)=T(n,1)+h/2*double(subs(f,x,a+(i-1/2)*h));
end
for i=2:n
k=4^(i-1);
T(n,i)=k/(k-1)*T(n,i-1)-1/(k-1)*T(n-1,i-1);
end
d=abs(T(n,n)-T(n-1,n-1));
end
I=T(n,n) %输出计算值
望采纳!谢谢!
8. 已知椭圆的方程是x^2+y^2/4=1,要在matlab用积分的方法求其周长,应该怎么做啊急啊,谢谢了
曲线长度积分l=∫√[φ'(t)^2+ψ'(t)^2] dt 椭圆参数方程为x=φ(t)=a sint y=ψ(t)=b cost
或l=∫√(1+y‘^2)dx 在(0,pi/2) 内 椭圆y=b*(1-x^2/a^2)^0.5
sin²x=(sinx)^2
9. 用matlab,龙贝格算法计算∫(0到1)[x/(4+x²)]dx的近似值。 求程序代码!!!
%龙贝格求积算法function I=romberg(a,b)h=b-a;T(1)=h/2*(fun(a)+fun(b));m=1;while 1 h=h/2; S(1)=1/2*T(1)+h*sumf(2^(m-1),a,h); for j=1:m S(j+1)=S(j)+(S(j)-T(j))/(4^j-1); end if abs(S(m+1)-T(m))<1e-6 break; end T=S;m=m+1;endI=S(m+1);end
function f=sumf(m,a,h)for j=1:m y(j)=fun(a+(2*j-1)*h);endf=sum(y);end
function f=fun(x)f=x/(4+x^2);end结果:
>> I=romberg(0,1)
I =
0.111571775646293
>>