can anyone help me run this, or get this running

Question

Relly Syam 2021-8-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/890332-can-anyone-help-me-run-this-or-get-this-running

回答： Walter Roberson 2021-8-2

clear; clc; 
format long e 
tic %bagian 
n=10 ;
syms c0 c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 t r0 r1 r2 r3 r4 r5 r6 r7 r8 r9 r10 r11 c0=0; 
c1=1/11; c2=2/11; c3=3/11; c4=4/11; c5=5/11; c6=6/11; c7=7/11; c8=8/11; c9=9/11; c10=10/11; c11=1; 
r0=0; r1=1; r2=2; r3=3; r4=4; r5=5; r6=6; r7=7; r8=8; r9=9; r10=10; r11=11; 
EvalAt = [c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11]; 
ktemp = arrayfun(@(EA) euler([r0, r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11], EA).', EvalAt, 'uniform', 0); 
ptemp=arrayfun(@(EA) int((3^(EA-t))*euler([r0, r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11],t),t,0,EA).', EvalAt, 'uniform', 0); 
E = horzcat(ktemp{:}).' ; 
K = horzcat(ptemp{:}).' ; 
Ek=E-K ; 
Inv_Ek= inv(Ek) ; 
F=[c0*3^c0;c1*3^c1;c2*3^c2;c3*3^c3;c4*3^c4;c5*3^c5;c6*3^c6;c7*3^c7;c8*3^c8;c9*3^c9;c10*3^c10;c11*3^c11]; 
C=Ek\F ; 
Ua=@(x)(C(1)*euler(0,x)+C(2)*euler(1,x)+C(3)*euler(2,x)+C(4)*euler(3,x)+C(5)*euler(4,x)+C(6)*euler(5,x)+C(7)*euler(6,x)+C(8)*euler(7,x)+C(9)*euler(8,x)+C(10)*euler(9,x)+C(11)*euler(10,x)+C(12)*euler(11,x)) ; 
Ue=@(x)((3^EA)*(1-e^(-EA))); 
uaa=zeros(12,1) ; 
uee=zeros(12,1) ; 
xx=zeros(12,1) ; 
k=0; 
for i=1:12 
    uaa(i)=Ua(k); 
    uee(i)=Ue(k); 
    xx(i)=k; 
    k=k+.1; 
end 
Uap= uaa; 
Uex= uee; 
y=(abs(uaa-uee)); 
[xx uee uaa y]; 
uee uaa y
%gambar plot(xx,uee,'k',xx,uaa,':o') grid on legend({'eksak','aproksimasi N=10'},'Location','Northwest')
%plot(xx,y,'k') %legend({'error N=10'},'Location','Northwest') %grid on xlabel('x'),ylabel('Eksak, Aproksimasi'); toc

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2021-8-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/890332-can-anyone-help-me-run-this-or-get-this-running#answer_758752

在 MATLAB Online 中打开

format long e 
n=10 ;
syms c0 c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 t r0 r1 r2 r3 r4 r5 r6 r7 r8 r9 r10 r11
c0=0; 
c1=1/11; c2=2/11; c3=3/11; c4=4/11; c5=5/11; c6=6/11; c7=7/11; c8=8/11; c9=9/11; c10=10/11; c11=1; 
r0=0; r1=1; r2=2; r3=3; r4=4; r5=5; r6=6; r7=7; r8=8; r9=9; r10=10; r11=11; 
EvalAt = [c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11]; 
ktemp = arrayfun(@(EA) euler([r0, r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11], EA).', EvalAt, 'uniform', 0); 
ptemp=arrayfun(@(EA) int((3^(EA-t))*euler([r0, r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11],t),t,0,EA).', EvalAt, 'uniform', 0); 
E = horzcat(ktemp{:}).' ; 
K = horzcat(ptemp{:}).' ; 
Ek=E-K ;
%do not calculate inv(Ek) !!
%{
Inv_Ek= inv(Ek) ; 
%}
F=[c0*3^c0;c1*3^c1;c2*3^c2;c3*3^c3;c4*3^c4;c5*3^c5;c6*3^c6;c7*3^c7;c8*3^c8;c9*3^c9;c10*3^c10;c11*3^c11]; 
C = double(Ek)\F ;                  %<=== important change!!
Ua=@(x)(C(1)*euler(0,x)+C(2)*euler(1,x)+C(3)*euler(2,x)+C(4)*euler(3,x)+C(5)*euler(4,x)+C(6)*euler(5,x)+C(7)*euler(6,x)+C(8)*euler(7,x)+C(9)*euler(8,x)+C(10)*euler(9,x)+C(11)*euler(10,x)+C(12)*euler(11,x)) ; 
Ue = @(EA)((3^EA)*(1-exp(-EA)));    %<=== recheck this !!!
uaa=zeros(12,1) ; 
uee=zeros(12,1) ; 
xx=zeros(12,1) ; 
k=0; 
for i=1:12 
    uaa(i)=Ua(k); 
    uee(i)=Ue(k); 
    xx(i)=k; 
    k=k+.1; 
end 
Uap= uaa; 
Uex= uee; 
y=(abs(uaa-uee)); 
[xx, uee, uaa, y]; 
uee, uaa, y
uee = 12×1
                         0
     1.062131630309665e-01
     2.258127092915629e-01
     3.603635391073476e-01
     5.116123773682127e-01
     6.815088885983266e-01
     8.722292439851662e-01
     1.086202425097018e+00
     1.326139582081997e+00
     1.595066801044776e+00
uaa = 12×1
    -6.938893903907228e-17
     1.173836988981493e-01
     2.758082659616968e-01
     4.864398971975612e-01
     7.632359801833405e-01
     1.123618200810357e+00
     1.589305303725420e+00
     2.187343074410759e+00
     2.951377916832975e+00
     3.923231266365520e+00
y = 12×1
     6.938893903907228e-17
     1.117053586718289e-02
     4.999555667013389e-02
     1.260763580902136e-01
     2.516236028151277e-01
     4.421093122120309e-01
     7.170760597402537e-01
     1.101140649313741e+00
     1.625238334750978e+00
     2.328164465320744e+00
%gambar plot(xx,uee,'k',xx,uaa,':o') grid on legend({'eksak','aproksimasi N=10'},'Location','Northwest')
%plot(xx,y,'k') %legend({'error N=10'},'Location','Northwest') %grid on xlabel('x'),ylabel('Eksak, Aproksimasi'); toc