In[]:=
$Version
Out[]=
13.0.0 for Microsoft Windows (64-bit) (December 3, 2021)
(*Analyticalmethod-----------------------------------------------------------------------------------------*)
In[]:=
ClearAll["`*"];Remove["`*"];
In[]:=
SOL1=DSolve[y[x]==Integrate[(Sin[x]+Cos[t])*y[t],{t,0,1}]+Cos[x],y[x],x]//FullSimplify
Out[]=
y[x]
2(Cos[1]+Sin[1])+Cos[x](-1+4Cos[1]+Cos[2]-2Sin[2])+(2Cos[2]+4Sin[1]+Sin[2])Sin[x]
-1+4Cos[1]+Cos[2]-2Sin[2]
In[]:=
Plot[y[x]/.SOL1,{x,0,1}]
Out[]=
(*Numericmethod1-----------------------------------------------------------------------------------------*)
In[]:=
n=5;(*Increasenifyouneedmoreaccuracy*)a1=0;b1=1;B[x_,i_]:=Binomial[n,i]*((x-a1)^i*(b1-x)^(n-i))/(b1-a1)^nk[t_,x_]:=Sin[x]+Cos[t];a[x_]:=1;f[x_]:=Cos[x];λ=1;nIntegrate[most__]:=NIntegrate[most,Method{Automatic,SymbolicProcessing0}]K2[j_?NumericQ]:=nIntegrate[B[x,j]f[x],{x,a1,b1}]//Quiet;mid[x_?NumericQ,i_?NumericQ]:=nIntegrate[λk[t,x]B[t,i],{t,a1,b1}]//Quiet;K3[j_?NumericQ,i_?NumericQ]:=nIntegrate[(a[x]B[x,i]-mid[x,i])B[x,j],{x,a1,b1}]//Quiet;SOL=Table[Total@Table[M[i]K3[j,i],{i,0,n}]K2[j],{j,0,n}];func=Total@Table[M[i]B[x,i],{i,0,n}]/.NSolve[SOL,Table[M[i],{i,0,n}],Reals][[1]]
Out[]=
-1.57426-11.0783x-29.071-35.9502-21.0887(1-x)-4.73261
5
(1-x)
4
(1-x)
3
(1-x)
2
x
2
(1-x)
3
x
4
x
5
x
In[]:=
Plot[{y[x]/.SOL1,func},{x,0,1},PlotStyle{Red,{Dashed,Black}},PlotLegends{"Analytical","Numeric"}]
Out[]=
In[]:=
error=Abs;(*Absoluteerror*)
(y[x]/.SOL1)-func
(y[x]/.SOL1)
In[]:=
Plot[error,{x,0,1},PlotRangeAll,PlotLegends{"Absolute error"}]
Out[]=
In[]:=
(*Numericmethod2-----------------------------------------------------------------------------------------*)
In[]:=
nN=5;(*Increasenifyouneedmoreaccuracy*)xcol=Table[-1+2j/nN,{j,0,nN}];v[x_]:=Table[BernoulliB[n,x],{n,0,nN}];A=Array[B,{nN+1}];u[x_]=A.v[x];eqs=Table[u[xcol[[i]]]-Cos[xcol[[i]]]-A.NIntegrate[v[t]*(Sin[xcol[[i]]]+Cos[t]),{t,0,1},AccuracyGoal8],{i,Length[xcol]}];sol10=FindRoot[Table[eqs[[i]]0,{i,Length[eqs]}],Table[{B[i],1/10},{i,Length[A]}]];func2=u[x]/.sol10
Out[]=
-3.20703-3.1583-+x+0.316832-x++0.527607-+-0.0248598-+-2+-0.025848-+-+
1
2
1
6
2
x
x
2
3
2
x
2
3
x
1
30
2
x
3
x
4
x
x
6
5
3
x
3
5
4
x
2
5
x
In[]:=
Plot[{func,func2},{x,0,1},PlotStyle{Red,{Dashed,Black}},PlotLegends{"Analytical","Numeric"}]
Out[]=
In[]:=
errorB=Abs;(*Absoluteerror*)
func2-func
func2
In[]:=
Plot[errorB,{x,0,1},PlotRangeAll,PlotLegends{"Absolute error"}]
Out[]=
(*Numericmethod3quasi-Simpson'srule-----------------------------------------------------------------------------------------*)
(*From:"https://mathematica.stackexchange.com/questions/23388/formatting-results-of-integral-approximation-methods-comparison"*)
In[]:=
simpson[func_,{x_,a_,b_},n_]:=Module[{f,h=(b-a)/n,wt},wt[i_]:=3.+(-1)^(i-1);f[t_]:=func/.x->t;If[EvenQ[n],(h/3)*(f[a]+Sum[wt[i]*f[a+i*h],{i,1,n-1}]+f[b]),"n must be even"]];
In[]:=
nN=10;(*Increasenifyouneedmoreaccuracy,nmustbeeven*)xcol=Table[-1+2j/nN,{j,0,nN}];v[x_]:=Table[BernoulliB[n,x],{n,0,nN}];A=Array[B,{nN+1}];u[x_]=A.v[x];eqs=Table[u[xcol[[i]]]-Cos[xcol[[i]]]-A.simpson[v[t]*(Sin[xcol[[i]]]+Cos[t]),{t,0,1},nN],{i,Length[xcol]}];sol11=FindRoot[Table[eqs[[i]]0,{i,Length[eqs]}],Table[{B[i],1/10},{i,Length[A]}]];func3=u[x]/.sol11
Out[]=
-3.20703-3.15832-+x+0.316397-x++0.526387-+-0.026367-+-2+-0.0263211-+-++0.000875914-+-3++0.000623297-+-+-0.0000182487-+-+-4+-0.0000100177-+2-+6-+-2.71017×-+5-7+-5+
1
2
1
6
2
x
x
2
3
2
x
2
3
x
1
30
2
x
3
x
4
x
x
6
5
3
x
3
5
4
x
2
5
x
1
42
2
x
2
5
4
x
2
5
x
6
x
x
6
7
3
x
6
7
5
x
2
7
6
x
2
7
x
1
30
2
2
x
3
7
4
x
3
14
6
x
3
7
x
8
x
3x
10
3
x
21
5
x
5
7
x
9
8
x
2
9
x
-7
10
5
66
3
2
x
2
4
x
6
x
15
8
x
2
9
x
10
x
In[]:=
Plot[{func,func3},{x,0,1},PlotStyle{Red,{Dashed,Black}},PlotLegends{"Analytical","Numeric"}]
Out[]=