试用matlab求如下数值积分(其中R,r为常量)
syms x y k kp delte
y = exp(-(k-kp)^2/σ^2)*cos(kx);
f = int(y, k, kp-2*delte, kp+2*delte);
MATLAB中主要用int进行符号积分,用trapz,dblquad,quad,quad8等进行数值积分。
int(s) 符号表达式s的不定积分
int(s,x) 符号表达式s关于变量x的不定积分
int(s,a,b) 符号表达式s的定积分,a,b分别为积分的上、下限
int(s,x,a,b) 符号表达式s关于变量x的定积分,a,b分别为积分的上、下限
trapz(x,y) 梯形积分法,x时表示积分区间的离散化向量,y是与x同维数的向量,表示被积函数,z返回积分值。
quad8(‘fun’,a,b,tol) 变步长数值积分,fun表示被积函数的M函数名,a,b分别为积分上、下限,tol为精度,缺省至为1e-3.
fblquad(‘fun’,a,b,c,d) 矩形区域二重数值积分,fun表示被积函数的M函数名,a,b分别为x的上、下限,c,d分别为y的上、下限.
例1 计算二重积分
先编写四个M函数文件,
%二重积分算法文件dblquad2.m
function S=dblquad2(f_name,a,b,c_lo,d_hi,m,n)
%其中f_name为被积函数字符串,'c_lo'和'd_hi'是y的下限和上限函数 ,都是x的标量函数;a,b分别为x的下限和上限;m,n分别为x和y方向上的等分数(缺省值为100).
if nargin<7, n=100; end
if nargin<6, m=100; end
if m<2|n<2
error('Numner of intervals invalid');
end
mpt=m+1; hx=(b-a)/m; x=a+(0:m)*hx;
for i=1:mpt
ylo=feval_r(c_lo,x(i)); yhi=feval_r(d_hi,x(i));
hy=(yhi-ylo)/n;
for k=1:n+1 y(i,k)=ylo+(k-1)*hy; f(i,k)=feval_r(f_name,x(i),y(i,k)); end
G(i)=trapz(y(i,:),f(i,:));
end
S=trapz(x,G);
%被积函数eg3_fun.m
function z=eg3_fun(x,y)
z=1+x+y;
%积分下限函数eg3_low.m
function y=eg3_low(x)
y=-sqrt(1-x^2);
%积分上限函数eg3_up.m
function y=eg3_up(x)
y=sqrt(1-x^2);
保存后,在命令窗口用MATLAB代码:
>>clear;
>>dblquad2('eg3_fun',-1,1,'eg3_low','eg3_up')
结果为
ans =3.1383
为了得到更精确的数值解,需将区间更细化,比如x和y方向等分为1000分,MATLAB代码:
>>clear; dblquad2('eg3_fun',-1,1,'eg3_low','eg3_up',1000,1000)
结果为 ans =3.1415。
此题也可用int符号计算求解,MATLAB代码为:
>>clear; syms x y;
>>iy=int(1+x+y,y,-sqrt(1-x^2),sqrt(1-x^2));
>>int(iy,x,-1,1)
结果为
ans =pi
例2 quad8计算定积分
%M函数fun1.m
function y=fun1(x)
y=x.^4;
保存后,在命令窗口用MATLAB代码:
>>clear;
>>quad8('fun1',-2,2)
>>vpa(quad8('fun1',-2,2),10) %以10位有效数字显示结果
结果为
ans =12.8000
ans =12.80000000
对于变步长数值积分,常用的有quad,quad8两种命令,quad使用自适应步长Simpson法, quad8使用自适应步长8阶Newton-Cotes法,我们建议用quad8,它不但精度较高,且对假收敛和假奇异积分具有一定的适应性,而quad较差..
龙贝格积分法MATLAB程序代码
function [I,step]=Roberg(f,a,b,eps)
if(nargin==3)
eps=1.0e-4;
end;
M=1;
tol=10;
k=0;
T=zeros(1,1);
h=b-a;
T(1,1)=(h/2)*(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b));
while tol>eps
k=k+1;
h=h/2;
Q=0;
for i=1:M
x=a+h*(2*i-1);
Q=Q+subs(sym(f),findsym(sym(f)),x);
end
T(k+1,1)=T(k,1)/2+h*Q;
M=2*M;
for j=1:k
T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1);
end
tol=abs(T(k+1,j+1)-T(k,j));
end
I=T(k+1,k+1);
step=k;
自适应法求积分MATLAB程序代码
function I=SmartSimpson(f,a,b,eps)
if(nargin==3)
eps=1.0e-4;
end;
e=5*eps;
I=SubSmartSimpson(f,a,b,e);
function q=SubSmartSimpson(f,a,b,eps)
QA=IntSimpson(f,a,b,1,eps);
QLeft=IntSimpson(f,a,(a+b)/2,1,eps);
QRight=IntSimpson(f,(a+b)/2,b,1,eps);
if(abs(QLeft+QRight-QA)<=eps)
q=QA;
else
q=SubSmartSimpson(f,a,(a+b)/2,eps)+SubSmartSimpson(f,(a+b)/2,b,eps);
end
线性振动响应分析的wilson θ积分法MATLAB代码
% see also http://www.matlabsky.com
% contact me matlabsky@gmail.com
% 2010-02-26 16:52:12
%
clc
clear
% 结构运动方程参数
M=1500000;
K=3700000;
C=470000;
% 威尔逊参数θ
theta=1.4;
dt=0.02; % 时间间隔
tau=dt*theta;
% 数据处理
eqd=load('acc_ElCentro_0.34g_0.02s.txt'); % 加速激励,第一列是时间,第二列是加速度
n=size(eqd,1);
t=0:dt:(n-1)*dt;
xg=eqd(:,2)*9.8; % 对加速度进行处理
dxg=diff(xg)*theta; %
F=-M*xg;
% D2x 加速度; Dx 速度; x 位移
D2x=zeros(n,1);
Dx=zeros(n,1);
x=zeros(n,1);
for i=1:n-1
K_ba=K+3/tau*C+6/tau^2*M;
dF_ba=-M*dxg(i)+(M*6/tau+3*C)*Dx(i)+(3*M+tau/2*C)*D2x(i);
dx=dF_ba/K_ba;
dD2x=(dx*6/tau^2-Dx(i)*6/tau-3*D2x(i))/theta;
D2x(i+1)=D2x(i)+dD2x;
Dx(i+1)=Dx(i)+D2x(i)*dt+dD2x/2*dt;
x(i+1)=x(i)+Dx(i)*dt+D2x(i)*dt^2/2+dD2x/6*dt^2;
end
subplot(311)
plot(t,x) % 位移
subplot(312)
plot(t,Dx) % 速度
subplot(313)
plot(t,D2x)% 加速度
常微分方程求解方法之四阶龙格-库塔算法matlab程序代码
function [x,y] = MyRunge_Kutta(fun,x0,xt,y0,PointNum,varargin)
%Runge-Kutta 方法解微分方程形为 y’(t) = f(x,y(x))
%此程序可解高阶的微分方程。只要将其形式写为上述微分方程的向量形式
% x范围为[x0,xt],初值为 y0, PointNum为离散点数,varargin为可选输入项可传适当参数给函数f(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]‘*h; %得x向量值
for k = 1:PointNum %迭代计算
f1 = h*feval_r(fun,x(k),y(k,:),varargin {:});
f1 = f1(:)’; %得公式中k1
f2 = h*feval_r(fun,x(k) + h/2,y(k,:) + f1/2,varargin{:});
f2 = f2(:)’; %得公式中k2
f3 = h*feval_r(fun,x(k) + h/2,y(k,:) + f2/2,varargin{:});
f3 = f3(:)’; %得公式中k3
f4 = h*feval_r(fun,x(k) + h,y(k,:) + f3,varargin{:});
f4 = f4(:)’; %得公式中k4
y(k + 1,:) = y(k,:) + (f1 + 2*(f2 + f3) + f4)/6; %
题主的做法主要存在两点问题:
1、quad函数用于计算数值积分,函数表达式中不能包含符号量;
2、被积函数的表达式应该写成关于被积变量的向量化的形式(也就是应该用点运算)。
参考代码:
R=1;
syms L;
rr = 0 : 0.1 : 1;
for ii = 1 : length(rr)
r = rr(ii);
f = @(l)(acos((1+l*l-r*r)/(2*l))+r*r*acos((r*r+l*l-1)/(2*r*l))-0.5*sqrt(4*r*r-(1+r*r-l*l)^2))*2*l/(pi*r^4);
fun = @(L) arrayfun(f,L);
J(ii) = quadl(fun,0,r);
end
plot(rr, J)
或者也可以借用楼上 枫箫1 的部分代码,写成:
R=1;
syms L;
rr = 0 : 0.1 : 1;
for ii = 1 : length(rr)
r = rr(ii);
SOA=R^2*acos((R^2+L^2-r^2)/(2*R*L))+r^2*acos((r^2+L^2-R^2)/(2*r*L))-...
0.5*sqrt(4*R^2*r^2-(R^2+r^2-L^2)^2);
PAB=SOA/(pi*r^2);
p=2*L/r^2;
f=PAB*p;
fun = eval(['@(L)' vectorize(f)]);
fun = @(l) arrayfun(@(L)eval(f),l);
J(ii) = quadl(fun,0,r);
end
plot(rr, J)
以上代码虽可以运行,但被积函数存在问题——SOA的第一项反余弦的值可能出现复数(因为在r稍小的情况下,acos的参数大于1),请题主再仔细检查一下。
R=10;
r=5;
syms L;
SOA=R^2*acos((R^2+L^2-r^2)/(2*R*L))+r^2*acos((r^2+L^2-R^2)/(2*r*L))-...
0.5*sqrt(4*R^2*r^2-(R^2+r^2-L^2)^2);
PAB=SOA/(pi*r^2);
p=2*L/r^2;
f=PAB*p;
fun=@(L) eval(f);
quadl(fun,0,r)
常量有具体的取值?
如何使用matlab,求AX=B?
MATLAB是matrix&laboratory两个词的组合,意为矩阵工厂。是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中,为科学研究、工程设计以及必须进行有效数值...
怎样用matlab求出sinx的值?
运行了你的程序,存在着下列几个问题:1、用下列这个代码求解,其t、y值是无法代入方程的。x=solve('465*sin(x)+200=120*cos(t*1.2+pi\/3)+369*sin(x-y)','x');所以软件会提示有关sym的信息。2、如要显示x值,其该语句最后不能有分号(;)。由于你给出的t、y值是以一组向量形式...
如何利用matlab求相关系数?
1、第一步我们首先需要知道matlab中求相关系数用到的是corrcoef函数,在命令行窗口中输入“help corrcoef”,可以看到corrcoef函数用法,2、第二步在命令行窗口中输入a=[1 3 6 7 8 16],b=[2 4 7 9 15 19],创建两个矩阵,求两个矩阵的相关系数,3、第三步输入corrcoef(a,b),按回车键,...
如何用matlab求解热传导偏微分方程?
1、首先要打开MATLAB R2016a软件,如下图所示。2、然后在打开的页面中,选择默认模式(Generic Scalar)-标量模式,具体如图。3、建立几何模型,绘制两个椭圆,再定义边界条件,具体如图所示。4、再定义PDE类型和系数,如下图所示。5、并将其三角形网格化,具体如图所示。6、最后可以对PDE图形进行求解了...
matlab如何求方差
matlab的方差求算在matlab程序上输入下列例子:Matlab 函数:var >>X=[1,2,3,4]>>var(X)=1.6667 >> sum((X(1,:)-mean(X)).^2)\/length(X)=1.2500 >> sum((X(1,:)-mean(X)).^2)\/(length(X)-1)=1.6667 var没有求矩阵的方差功能,可使用std先求均方差,再平方得到方差。st...
如何用matlab求解一个二阶常系数微分方程组
2、 输入微分方程求解程序-->点击保存-->点击运行。3、在matlab的命令窗口即可看到求解结果,是一个关于参数a,b的表达式 第二种方法:利用Matlab中的solver函数(包括ode45、ode23、ode15s等)来求解微分方程的数值解,这种方法是最常用的方法,对于dsolve函数难以求解的方程就可以利用这种方法求解方程的...
如何用MATLAB求逆矩阵
求P,Q的交集,这一步有专门的凸集分离定理Farkas定理。如何用matlab 求矩阵的逆 可以调用matlab中的 inv 函数。调用格式如下:Y=inv(x)输入矩阵X必须为方阵。输出Y的精度默认为0.0001.如何用cublas计算逆矩阵 一般考试的时候,矩阵求逆最简单的办法是用增广矩阵 如果要求逆的矩阵是A 则对增广...
求用Matlab如何画求导函数曲线
具体如下:1、第一步,打开matlab软件,出现如下界面,见下图,转到下面的步骤。2、第二步,完成上述步骤后,敲入命令“clear;clc; ”来清理工作空间,见下图,转到下面的步骤。3、第三步,完成上述步骤后,敲入命令“syms x”来定义一个符号变量,见下图,转到下面的步骤。4、第四步,完成上述...
如何用matlab求解点如: A 418 45 B 570 136 C 477 265 D 417 194 E...
>> X = [418 570 477 417 338 206 418];>> Y = [45 136 265 194 239 173 45];>> S = polyarea(X, Y)S = 39848>> plot(X, Y);
MATLAB中如何求幂函数,比如0.9的n次方,n=[-5:5]
1、首先双击matlab软件图标,打开matlab软件,可以看到matlab软件的界面。2、使用“0.1:0.1:5;”创建一维数组,表示从0.1到5,每隔0.1会取一个数字。这个一维数组用来作为一元一次函数的横坐标的数值。3、接着创建三个幂函数,分别是y1=x.^(1\/4); y2=x.^(1\/2); y3=x.^(3\/2)。4、...
丙疯天君: 题主的做法主要存在两点问题:1、quad函数用于计算数值积分,函数表达式中不能包含符号量;2、被积函数的表达式应该写成关于被积变量的向量化的形式(也就是应该用点运算). 参考代码:R=1; syms L; rr = 0 : 0.1 : 1; for ii = 1 : length(rr) ...
都匀市13486617161: 如何用Matlab求以下积分 - ?
丙疯天君: 求数值积分可用quad函数.function test() clear all;clc; global alpha1 alpha2 m1 n1 m2 n2 E01 E02 b1 b2 R1 R2 alpha1=1; alpha2=1; m1=.2; n1=.3; m2=.4; n2=.5; E01=.6; E02=.7; b1=.8; b2=.9; R1=1; R2=1; u=0; v=.01; Q = quad(@model,u,v) end ...
都匀市13486617161: 如何用matlab 数值法算这个积分 - ?
丙疯天君: 程序: fun=sin(0.5*pi*x./y);%% a=int(int(fun,y,sqrt(x),x),x,1,2); b=simple(a) %化简 I=vpa(b,4) %得到4位近似解,也可以任意N位解 结果: I = 0.2719
都匀市13486617161: 用MATLAB算积分 - ?
丙疯天君: 使用MATLAB软件,可以用int()计算不定积分或定积分.计算方法如下:syms x int((log(x))^2/x) %这里 ln(x) 用 log(x) 来表示 计算结果 log(x)^3/3 %ln³(x)/3
都匀市13486617161: matlab用quad函数求数值积分 - ?
丙疯天君: 应这样改: 一个文件中内容如下: function y = fun(x) y = x.^2+3*x.^3-1; end 另一个主文件内容如下: quad('fun',-1,1)
都匀市13486617161: matlab 求解如下积分 - ?
丙疯天君: 使用数值解法吧 x=1%给x赋值 fun=@(w)exp(i*x*cos(w)) quadl(fun,0,2*pi)
都匀市13486617161: matlab 求解积分函数∫(0,1)e²x dx求解函数 的数值积分和符号积分并比较结果求解积分函数∫(0,1)e²x dx求解函数 的数值积分和符号积分并比较 - ?
丙疯天君:[答案] 符号 syms x; int(exp(2*x),x,0,1) ans = exp(2)/2 - 1/2 数值 f=@(x)exp(2*x);quad(f,0,1) ans = 3.1945 符号积分精确度高但速度慢,有时候有些函数没有解析解,就得用数值积分,并且数值积分速度快,但精确度不高
都匀市13486617161: 大家帮忙用matlab进行数值积分,并列出命令行 - ?
丙疯天君: 用matlab可以如下数值积分法,来求解定积分、二重积分、三重积分的数值解问题. 1、梯形数值积分计算 trapz() X = 0:pi/100:pi; Y = sin(X); Z = pi/100*trapz(Y) 2、自适应辛普森数值积分计算 quad() F = @(x)1./(x.^3-2*x-5); Q = quad(F,0,2); 3、自...
都匀市13486617161: matlab中怎样对二元函数中的一个变量做数值积分? - ?
丙疯天君: clear all clc syms va rho tau f=@(va,rho,tau) (2.5-0.1*(va+rho*sin(tau)))./(exp(2.5-0.1*(va+rho*sin(tau)))-1); n = 5000;%当想要结果更精确时,可以把n设置更大 tau = linspace(0,2*pi,n);%自编的简单方法 Tn =@(va,rho) pi/n*(sum(f(va,rho,tau(1:end-...
都匀市13486617161: matlab中: 求解f(x)=1/(1+x^i ) 分别求出 i 从1到20在区间[ - 1,1]上的数值积分,求解决办法? - ?
丙疯天君:[答案] 具体代码如下: for i=1:20 f=@(x)1./(1+x.^i); result(i)=quad(f,-1,1); end 不过我要说一点的是,当i=1的时候,积分就是无穷了哦~