matlab 符號積分和數值積分

2021-03-12 16:08:03 字數 5205 閱讀 6625

1樓:楊好巨蟹座

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');

endmpt=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,:));

ends=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);

endt(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);

endtol=abs(t(k+1,j+1)-t(k,j));

endi=t(k+1,k+1);

step=k;

自適應法求積分matlab程式**

function i=**artsimpson(f,a,b,eps)

if(nargin==3)

eps=1.0e-4;

end;

e=5*eps;

i=sub**artsimpson(f,a,b,e);

function q=sub**artsimpson(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=sub**artsimpson(f,a,(a+b)/2,eps)+sub**artsimpson(f,(a+b)/2,b,eps);

end線性振動響應分析的wilson θ積分法matlab**

% 結構運動方程引數

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;

endsubplot(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;

endif nargin < 3

y0 = 0;

endy(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; %

matlab求積分的數值解,MATLAB求積分的數值解

可以積分出來的,只是好像你的r和 怎麼是反的 clear all clc syms x r 0.3136 beta 0.579 y r 2 pi beta 0.5 x 1.5 exp x r 2 2 beta x s int x exp 3 x y,x,0,inf s double s s 0.08...

定積分中符號的問題,定積分中符號的問題

解 dx就是對x積分啊 在數學上表示 且 的意思 積分上限就是那個大的數,積分下限就是小的數,就是在小的數到大的數的區間內積分啊,這就是上限和下限。倒寫的v是次方的意思 沒槓的f上面是積分上限,下面表示積分下限 在一個線段上插入若干分點 x1小於x2小於 小於xn 1 小於xn 沒看懂什麼意思,最好...

這個積分符號是什麼意思,積分符號上加個圈是什麼意思

表示在閉曲線上的曲線積分。在閉曲線上的曲線積分可以寫成 也可以仍然寫 但不是閉曲線上的曲線積分,不可以用 只能寫 積分符號上加個圈是什麼意思?10 積分符號上加個圈表示強調。在閉曲線上的曲線積分,我們常常在積分號上面寫一個 以強調,即把積分號 寫成 在閉曲線上的曲線積分可以寫 也可以仍然寫 但不是閉...