Fredholm第⼆类积分⽅程的MATLAB代码实现(1)
这是《The Numerical Solution of Integral Equations of the Second Kind》书上的例题代码实现,基于MATLAB来编写的,对应于这⼀篇的Python代码
本代码主要实现Fredholm integral equations of the second kind,⽅程如下:
lamda*x(t)-integrate(K(t,s)*x(s))ds=y(t),a=0<=t<=b,K(t,s)取exp(s*t)
⾸先,我们给定两个精确解exp(-t)cos(t)和sqrt(t),求出其对应的y(t),然后再来反解x(t).更详细说明可
参见《The Numerical Solution of Integral Equations of the Second Kind》P30-34.
⾸先定义各种函数,保存为.m⽂件
第⼀步:定义核函数
function k=func_k(s,t)
k=exp(s*t); %k的表达式
第⼆步:定义给定的精确解x
function x=func_x(t)
x=exp(-t)*cos(t); %x的表达式
第三步:⽤精确解求出y
function y=func_y(x,k,a,b,lamda)
syms t s;
y=subs(lamda*x-int(x*k,t,a,b),t,s); % ⽤符号积分求出y的表达式,并统⼀y中的变量为s
第四步:求⽅程组的右端项(相当于Ax=b中的b)
function f=func_f(y,a,b,lamda,n)
matlab定义函数表达式syms s;
f=ones(n,1);
for i=1:n
f(i)=int(y*s^(i-1),s,a,b); %这⾥⽤符号积分求y*s^(i-1)
end
第五步:求系数c
function c=solve_c(f,a,b,lamda,n) %
B=ones(n,n); %⽤于存放系数矩阵
for i=1:n
for j=1:n
B(i,j)=-b^(i+j-1)/(factorial(j-1)*(i+j-1));
end
B(i,i)=lamda+B(i,i);
end
c=B\f; %求c,其为⼀个n维的向量
end
第六步:求xn
function xn=solve_xn(y,c,a,b,lamda,n)
syms t;
m=0;
for i=0:n-1
m=m+c(i+1)*t^i/factorial(i);
end
xn=(y+m)/lamda;
第七步:求⽆穷范数
function E=err(x,xn,a,b)
syms s t;
tt=a:0.1:b;
E=vpa(norm(subs(subs(x-xn,s,t),t,tt),Inf),20); %取区间[a,b]中的离散点,求⽆穷范数
以上是定义好的函数,接下来再写⼀个主函数⽂件,然后运⾏此⽂件即可。以上的核函数,x函数等可以改成⾃⼰想要的
clear
n=input('please input n:')
a=input('please input a:')
b=input('please input b:')
lamda=input('please input lamda:')
e=ones(n,1); %⽤来存误差
syms s t;
for i=1:n
k=func_k(s,t);
x=func_x(t);
y=func_y(x,k,a,b,lamda);
f=func_f(y,a,b,lamda,i);
c=solve_c(f,a,b,lamda,i);
xn=solve_xn(y,c,a,b,lamda,i);
e(i)=err(x,xn,a,b);
end
e
if n==10
figure(1)
tt=a:0.1:b;
plot(tt,vpa(subs(subs(x-xn,s,t),t,tt),20),'r--')
xlabel('t'), ylabel('error')
title('error,n=10,a=0,b=1,lamda=5')
end
运⾏结果我就不放上来了,因为这台电脑没有按MATLAB,o(╥﹏╥)o
版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。
发表评论