From 67448b9e87104c681ce6149ea1dda5c0eb67e601 Mon Sep 17 00:00:00 2001 From: geniusjoe <843520313@qq.com> Date: Thu, 14 Mar 2019 23:48:59 +0800 Subject: [PATCH 01/19] Create newton_iteration.m --- newton_iteration.m | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100644 newton_iteration.m diff --git a/newton_iteration.m b/newton_iteration.m new file mode 100644 index 0000000..81ca7c3 --- /dev/null +++ b/newton_iteration.m @@ -0,0 +1,28 @@ +function res=newton_iteration(first_x_value,fx) +% \param x: start position +% \param fx: syms expression +% \return res: the value the root in +% +% can change precision and max_iteration_times at the following +disp(fx); +precision=1e-6; +max_iteration_times=100000; +% ------------------------------------------------------------% +syms x; +diff_fx=diff(fx); +f_x_next=(x-fx/diff_fx); + +x_before=first_x_value; +x_next=subs(f_x_next,x,first_x_value); + +iterate_times=1; +while (abs(x_next-x_before)>precision) && (iterate_times<=max_iteration_times) + current_display=[x_next;iterate_times]; + formatSpec = '\nthe current x_value is %06f,\nthe current iterate_times is %d'; + fprintf(formatSpec,current_display); + x_before=x_next; + x_next=double(subs(f_x_next,x,x_before)); + iterate_times=iterate_times+1; +end +res=x_next; +end \ No newline at end of file From 93c3e9bb98bc4381e0defef9f8092cd8fa2cc3b3 Mon Sep 17 00:00:00 2001 From: geniusjoe <843520313@qq.com> Date: Fri, 15 Mar 2019 00:39:12 +0800 Subject: [PATCH 02/19] Update newton_iteration.m --- newton_iteration.m | 27 +++++++++++++++++++-------- 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/newton_iteration.m b/newton_iteration.m index 81ca7c3..4ee4261 100644 --- a/newton_iteration.m +++ b/newton_iteration.m @@ -3,25 +3,36 @@ % \param fx: syms expression % \return res: the value the root in % -% can change precision and max_iteration_times at the following -disp(fx); +% can change precision and max_iteration_times and min_diff_value at the following precision=1e-6; -max_iteration_times=100000; +max_iteration_times=10; +min_diff_value=1e-4; % ------------------------------------------------------------% syms x; diff_fx=diff(fx); -f_x_next=(x-fx/diff_fx); +fx_next=x-fx/diff_fx; -x_before=first_x_value; -x_next=subs(f_x_next,x,first_x_value); +x_before=double(first_x_value); +x_next=double(subs(fx_next,x,x_before)); iterate_times=1; while (abs(x_next-x_before)>precision) && (iterate_times<=max_iteration_times) + % if f'x < min_diff_value, + % throw exception + if abs(subs(diff_fx,x,x_before)) Date: Sat, 16 Mar 2019 19:41:57 +0800 Subject: [PATCH 03/19] change disp() function into fprintf function --- newton_iteration.m | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/newton_iteration.m b/newton_iteration.m index 4ee4261..4aec4d7 100644 --- a/newton_iteration.m +++ b/newton_iteration.m @@ -5,7 +5,7 @@ % % can change precision and max_iteration_times and min_diff_value at the following precision=1e-6; -max_iteration_times=10; +max_iteration_times=20; min_diff_value=1e-4; % ------------------------------------------------------------% syms x; @@ -18,17 +18,16 @@ iterate_times=1; while (abs(x_next-x_before)>precision) && (iterate_times<=max_iteration_times) % if f'x < min_diff_value, - % throw exception + % throw warning if abs(subs(diff_fx,x,x_before)) Date: Sun, 17 Mar 2019 00:12:29 +0800 Subject: [PATCH 04/19] initial legendre --- Legendre.m | 33 +++++++++++++++++++++++++++++++++ newton_iteration.m | 20 +++++++++++++++----- 2 files changed, 48 insertions(+), 5 deletions(-) create mode 100644 Legendre.m diff --git a/Legendre.m b/Legendre.m new file mode 100644 index 0000000..7560834 --- /dev/null +++ b/Legendre.m @@ -0,0 +1,33 @@ +function res=Legendre() +% P1(X)=1 +% P2(X)=X +% ...... +% pN+2(X)= (2*n + 3) / (n + 2) * x * PN+1(x) - (n + 1) / (n + 2) * PN(x) +% +% can change precision and max_iteration_times and min_diff_value at the following +max_iteration_times=20; +% n=1...max_iteration_times +% ------------------------------------------------------------% +% example: +% >> syms x; +% fx=x-exp(-x); +% newton_iteration(0.5,fx); +% +% return +% the current x_value is 0.567143, +% the current iterate_times is 2 +% -- + +syms x; +coefficient_list=zeros(0,'syms',max_iteration_times); +coefficient_list(1)=1; +coefficient_list(2)=x; + +for n=1:max_iteration_times-2 + coefficient_list(n+2)= simplify(... + (2 * n + 3) / (n + 2) * x * coefficient_list(n+1) ... + - (n + 1) / (n + 2) * coefficient_list(n) ... + ); +end +pass=1 +end \ No newline at end of file diff --git a/newton_iteration.m b/newton_iteration.m index 4aec4d7..5619838 100644 --- a/newton_iteration.m +++ b/newton_iteration.m @@ -1,13 +1,23 @@ function res=newton_iteration(first_x_value,fx) % \param x: start position % \param fx: syms expression -% \return res: the value the root in -% +% \return res: the root value of fx +% % can change precision and max_iteration_times and min_diff_value at the following precision=1e-6; max_iteration_times=20; min_diff_value=1e-4; % ------------------------------------------------------------% +% example: +% >> syms x; +% fx=x-exp(-x); +% newton_iteration(0.5,fx); +% +% return +% the current x_value is 0.567143, +% the current iterate_times is 2 +% ------------------------------------------------------------% + syms x; diff_fx=diff(fx); fx_next=x-fx/diff_fx; @@ -20,14 +30,14 @@ % if f'x < min_diff_value, % throw warning if abs(subs(diff_fx,x,x_before)) Date: Sun, 17 Mar 2019 10:09:17 +0800 Subject: [PATCH 05/19] finish legendre.m --- Legendre.m | 34 +++++++++++++++++----------------- newton_iteration.m | 2 +- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/Legendre.m b/Legendre.m index 7560834..671ba76 100644 --- a/Legendre.m +++ b/Legendre.m @@ -1,33 +1,33 @@ -function res=Legendre() +function res=Legendre(first_x_value,iterate_times) +% \param first_x_value: start position +% \param iterate_times: max N num of PN +% \return res: the root value of Legendre +% ------------------------------------------------------------% % P1(X)=1 % P2(X)=X % ...... % pN+2(X)= (2*n + 3) / (n + 2) * x * PN+1(x) - (n + 1) / (n + 2) * PN(x) -% -% can change precision and max_iteration_times and min_diff_value at the following -max_iteration_times=20; -% n=1...max_iteration_times % ------------------------------------------------------------% % example: -% >> syms x; -% fx=x-exp(-x); -% newton_iteration(0.5,fx); +% >> Legendre(0.93,7); % -% return -% the current x_value is 0.567143, +% return +% the current x_value is 0.932470, % the current iterate_times is 2 -% -- +% ------------------------------------------------------------% syms x; -coefficient_list=zeros(0,'syms',max_iteration_times); +coefficient_list=sym(zeros(0,iterate_times)); coefficient_list(1)=1; coefficient_list(2)=x; -for n=1:max_iteration_times-2 - coefficient_list(n+2)= simplify(... - (2 * n + 3) / (n + 2) * x * coefficient_list(n+1) ... - - (n + 1) / (n + 2) * coefficient_list(n) ... +for index=1:iterate_times-2 + n=index-1; + coefficient_list(index+2)= simplify(... + (2 * n + 3) / (n + 2) * x * coefficient_list(index+1) ... + - (n + 1) / (n + 2) * coefficient_list(index) ... ); end -pass=1 +disp(transpose(coefficient_list)); +res=newton_iteration(first_x_value,coefficient_list(iterate_times)); end \ No newline at end of file diff --git a/newton_iteration.m b/newton_iteration.m index 5619838..9f55e38 100644 --- a/newton_iteration.m +++ b/newton_iteration.m @@ -1,5 +1,5 @@ function res=newton_iteration(first_x_value,fx) -% \param x: start position +% \param first_x_value: start position % \param fx: syms expression % \return res: the root value of fx % From 8fe9d895de32b33209826c8dbba810f6975b298a Mon Sep 17 00:00:00 2001 From: geniusjoe <843520313@qq.com> Date: Sun, 17 Mar 2019 13:20:59 +0800 Subject: [PATCH 06/19] Update newton_iteration.m --- newton_iteration.m | 67 +++++++++++++++++++++++++--------------------- 1 file changed, 37 insertions(+), 30 deletions(-) diff --git a/newton_iteration.m b/newton_iteration.m index 9f55e38..ab27939 100644 --- a/newton_iteration.m +++ b/newton_iteration.m @@ -1,5 +1,5 @@ -function res=newton_iteration(first_x_value,fx) -% \param first_x_value: start position +function res=newton_iteration(x_value_list,fx) +% \param x_value_list: list of start position % \param fx: syms expression % \return res: the root value of fx % @@ -9,40 +9,47 @@ min_diff_value=1e-4; % ------------------------------------------------------------% % example: -% >> syms x; -% fx=x-exp(-x); -% newton_iteration(0.5,fx); +% syms x; +% fx=x^2-2*x*exp(-x)+exp(-2*x); +% newton_iteration([0.5 0.6 0.7],fx); % % return -% the current x_value is 0.567143, -% the current iterate_times is 2 +% the answer is +% 0.5671 +% 0.5672 +% 0.5672 % ------------------------------------------------------------% syms x; diff_fx=diff(fx); fx_next=x-fx/diff_fx; - -x_before=double(first_x_value); -x_next=double(subs(fx_next,x,x_before)); - -iterate_times=1; -while (abs(x_next-x_before)>precision) && (iterate_times<=max_iteration_times) - % if f'x < min_diff_value, - % throw warning - if abs(subs(diff_fx,x,x_before))precision) && (iterate_times<=max_iteration_times) + % if f'x < min_diff_value, + % return + if abs(subs(diff_fx,x,x_before)) Date: Sun, 17 Mar 2019 13:42:26 +0800 Subject: [PATCH 07/19] modify fprintf --- Legendre.m | 16 +++++++++++----- newton_iteration.m | 4 ++-- 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/Legendre.m b/Legendre.m index 671ba76..da9daba 100644 --- a/Legendre.m +++ b/Legendre.m @@ -1,6 +1,7 @@ -function res=Legendre(first_x_value,iterate_times) +function res=Legendre(x_value_list,iterate_times) % \param first_x_value: start position % \param iterate_times: max N num of PN +% % \return res: the root value of Legendre % ------------------------------------------------------------% % P1(X)=1 @@ -9,11 +10,16 @@ % pN+2(X)= (2*n + 3) / (n + 2) * x * PN+1(x) - (n + 1) / (n + 2) * PN(x) % ------------------------------------------------------------% % example: -% >> Legendre(0.93,7); +% >> Legendre([0.93 -0.93 0.66 -0.66 0.23 -0.23],7); % % return -% the current x_value is 0.932470, -% the current iterate_times is 2 +% the answer is +% 0.932470 +% -0.932470 +% 0.661209 +% -0.661209 +% 0.238619 +% -0.238619 % ------------------------------------------------------------% syms x; @@ -29,5 +35,5 @@ ); end disp(transpose(coefficient_list)); -res=newton_iteration(first_x_value,coefficient_list(iterate_times)); +res=newton_iteration(x_value_list,coefficient_list(iterate_times)); end \ No newline at end of file diff --git a/newton_iteration.m b/newton_iteration.m index ab27939..e27acd0 100644 --- a/newton_iteration.m +++ b/newton_iteration.m @@ -40,7 +40,7 @@ end % print each iterate x_next value - fprintf('\nthe current x_value is %06f,\nthe current iterate_times is %d\n', ... + fprintf('\nthe current x_value is %f,\nthe current iterate_times is %d\n', ... x_next,iterate_times); % next iterate @@ -51,5 +51,5 @@ res(i)=x_next; end fprintf('\nthe answer is\n'); -disp(transpose(res)); +fprintf('\n%f\n',res); end From cb6874f84a62483927a72bb9478314a0c127d210 Mon Sep 17 00:00:00 2001 From: geniusjoe <843520313@qq.com> Date: Sun, 17 Mar 2019 22:06:51 +0800 Subject: [PATCH 08/19] Update README.md --- README.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index a7765e4..1cc7cbb 100644 --- a/README.md +++ b/README.md @@ -15,8 +15,9 @@ |算法名称 | 文件名 | |:--------:|:--------:| |共轭梯度法 |`CG_equ.m`| +|牛顿迭代法 |`newton_iteration.m`| ## 插值 | 算法名称 | | 文件名 | |:----------:|:-----------:|:---------------:| -|三次样条插值 |第一种边界条件
第二种边界条件
第三种边界条件|`Spline1_inter.m`
`Spline2_inter.m`
`Spline3_inter.m` | \ No newline at end of file +|三次样条插值 |第一种边界条件
第二种边界条件
第三种边界条件|`Spline1_inter.m`
`Spline2_inter.m`
`Spline3_inter.m` | From 17d2e1caca15116a196c282e48c106f68ccfdb83 Mon Sep 17 00:00:00 2001 From: geniusjoe <843520313@qq.com> Date: Tue, 9 Apr 2019 23:35:50 +0800 Subject: [PATCH 09/19] add romberg.m --- .gitignore | 1 + Romberg_new.m | 28 ++++++++++++++++++++++++++++ 2 files changed, 29 insertions(+) create mode 100644 Romberg_new.m diff --git a/.gitignore b/.gitignore index 3e3461a..942356c 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ *.asv +Romberg.m diff --git a/Romberg_new.m b/Romberg_new.m new file mode 100644 index 0000000..4fba55f --- /dev/null +++ b/Romberg_new.m @@ -0,0 +1,28 @@ +function res=Romberg_new(fx,a,b) +eps=1e-5; +buf=zeros(10,10); +i=1; +h=(b-a)/2; +buf(1,1)=(b-a)./2.*(fx(b)-fx(a)); +buf(1,2)=h/2*(fx(b)+fx(a)+2*fx(a+h*i)); +buf(2,1)=(4*buf(1,2)-buf(1,1))/(4-1); +if (abs((buf(i+1,1)-buf(i,1)))=eps) + i=i+1; + for m=2:i+1 + k=i-m; + n=2.^i; + h=(b-a)/2.^i; + buf(1,i+1)=h/2*(fx(a)+fx(b)); + for j=1:n-1 + buf(1,i+1)=buf(1,i+1)+h/2*(2*fx(a+j*h)); % get T_1_i+1 value + end + buf(m+1,k+1)=(4^m*buf(m,k+2)-buf(m,k+1))/(4^m-1); + end + end +end +disp(buf); +res=buf(a+1,b+1); +end \ No newline at end of file From bfc534d78e08b78477b1ad09f2850cacbecb23c5 Mon Sep 17 00:00:00 2001 From: geniusjoe <843520313@qq.com> Date: Wed, 10 Apr 2019 00:12:47 +0800 Subject: [PATCH 10/19] Update Romberg_new.m --- Romberg_new.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Romberg_new.m b/Romberg_new.m index 4fba55f..8162737 100644 --- a/Romberg_new.m +++ b/Romberg_new.m @@ -2,7 +2,7 @@ eps=1e-5; buf=zeros(10,10); i=1; -h=(b-a)/2; +h=(b-a); buf(1,1)=(b-a)./2.*(fx(b)-fx(a)); buf(1,2)=h/2*(fx(b)+fx(a)+2*fx(a+h*i)); buf(2,1)=(4*buf(1,2)-buf(1,1))/(4-1); From 89adc7f454b9cb8828aba53b712d180b8e058f28 Mon Sep 17 00:00:00 2001 From: geniusjoe <843520313@qq.com> Date: Wed, 10 Apr 2019 00:17:34 +0800 Subject: [PATCH 11/19] change_name --- Romberg_calculus.m | 55 ++++++++++++++++++++++++++++++++++++++++++++++ Romberg_new.m | 28 ----------------------- 2 files changed, 55 insertions(+), 28 deletions(-) create mode 100644 Romberg_calculus.m delete mode 100644 Romberg_new.m diff --git a/Romberg_calculus.m b/Romberg_calculus.m new file mode 100644 index 0000000..38e2e86 --- /dev/null +++ b/Romberg_calculus.m @@ -0,0 +1,55 @@ +function res=Romberg_calculus(fx,a,b) +% \param fx: syms expression +% \param a: left bound +% \param b: right bound +% \return res: the root value of fx +% +% can change precision and max_iteration_times and min_diff_value at the following +precision=1e-6; +max_iteration_times=20; +min_diff_value=1e-4; +% ------------------------------------------------------------% +% example: +% fx=@(x)x^2*exp(x); +% a=0; +% b=0.8; +% Romberg_new(fx,a,b); +% +% return +% the answer is +% 0.5697 +% 2.0890 +% 0.2099 +% 0.3196 +% 0.3153 +% 0.3148 +% 0.3146 +% 0.3146 +% ------------------------------------------------------------% +eps=1e-6; +buf=zeros(10,10); +i=1; +h=(b-a); +buf(1,1)=(b-a)./2.*(fx(b)-fx(a)); +buf(1,2)=h/2*(fx(b)+fx(a)+2*fx(a+h*i)); +buf(2,1)=(4*buf(1,2)-buf(1,1))/(4-1); +if (abs((buf(i+1,1)-buf(i,1)))=eps) + i=i+1; + for m=2:i+1 + k=i-m; + n=2.^i; + h=(b-a)/2.^i; + buf(1,i+1)=h/2*(fx(a)+fx(b)); + for j=1:n-1 + buf(1,i+1)=buf(1,i+1)+h/2*(2*fx(a+j*h)); % get T_1_i+1 value + end + buf(m+1,k+1)=(4^m*buf(m,k+2)-buf(m,k+1))/(4^m-1); + end + end +end +disp(buf); +res=buf(a+1,b+1); +end \ No newline at end of file diff --git a/Romberg_new.m b/Romberg_new.m deleted file mode 100644 index 8162737..0000000 --- a/Romberg_new.m +++ /dev/null @@ -1,28 +0,0 @@ -function res=Romberg_new(fx,a,b) -eps=1e-5; -buf=zeros(10,10); -i=1; -h=(b-a); -buf(1,1)=(b-a)./2.*(fx(b)-fx(a)); -buf(1,2)=h/2*(fx(b)+fx(a)+2*fx(a+h*i)); -buf(2,1)=(4*buf(1,2)-buf(1,1))/(4-1); -if (abs((buf(i+1,1)-buf(i,1)))=eps) - i=i+1; - for m=2:i+1 - k=i-m; - n=2.^i; - h=(b-a)/2.^i; - buf(1,i+1)=h/2*(fx(a)+fx(b)); - for j=1:n-1 - buf(1,i+1)=buf(1,i+1)+h/2*(2*fx(a+j*h)); % get T_1_i+1 value - end - buf(m+1,k+1)=(4^m*buf(m,k+2)-buf(m,k+1))/(4^m-1); - end - end -end -disp(buf); -res=buf(a+1,b+1); -end \ No newline at end of file From e04774af248976785c2c85f3bb480b9995388a7f Mon Sep 17 00:00:00 2001 From: geniusjoe <843520313@qq.com> Date: Wed, 10 Apr 2019 00:23:14 +0800 Subject: [PATCH 12/19] Update Romberg_calculus.m --- Romberg_calculus.m | 30 ++++++++++++++---------------- 1 file changed, 14 insertions(+), 16 deletions(-) diff --git a/Romberg_calculus.m b/Romberg_calculus.m index 38e2e86..c5ed148 100644 --- a/Romberg_calculus.m +++ b/Romberg_calculus.m @@ -2,12 +2,11 @@ % \param fx: syms expression % \param a: left bound % \param b: right bound -% \return res: the root value of fx +% \return res: the fx calculas % % can change precision and max_iteration_times and min_diff_value at the following -precision=1e-6; -max_iteration_times=20; -min_diff_value=1e-4; +eps=1e-6; +buf=zeros(15); % size % ------------------------------------------------------------% % example: % fx=@(x)x^2*exp(x); @@ -26,8 +25,6 @@ % 0.3146 % 0.3146 % ------------------------------------------------------------% -eps=1e-6; -buf=zeros(10,10); i=1; h=(b-a); buf(1,1)=(b-a)./2.*(fx(b)-fx(a)); @@ -38,18 +35,19 @@ else while (abs((buf(i+1,1)-buf(i,1)))>=eps) i=i+1; + n=2.^i; + h=(b-a)/(2^i); + buf(1,i+1)=h/2*(fx(a)+fx(b)); + for j=1:n-1 + buf(1,i+1)=buf(1,i+1)+h/2*(2*fx(a+j*h)); % get T_1_i+1 value + end + for m=2:i+1 - k=i-m; - n=2.^i; - h=(b-a)/2.^i; - buf(1,i+1)=h/2*(fx(a)+fx(b)); - for j=1:n-1 - buf(1,i+1)=buf(1,i+1)+h/2*(2*fx(a+j*h)); % get T_1_i+1 value - end - buf(m+1,k+1)=(4^m*buf(m,k+2)-buf(m,k+1))/(4^m-1); + k=i+2-m; + buf(m,k)=(4^m*buf(m-1,k+1)-buf(m-1,k))/(4^m-1); end end end -disp(buf); -res=buf(a+1,b+1); +disp(buf(:,1)); +res=buf(i,1); end \ No newline at end of file From 9a786cdec9635978a9e26fddee0fb188fc72caf7 Mon Sep 17 00:00:00 2001 From: geniusjoe <843520313@qq.com> Date: Wed, 10 Apr 2019 01:02:57 +0800 Subject: [PATCH 13/19] Create Runge_Kutta.m --- Runge_Kutta.m | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) create mode 100644 Runge_Kutta.m diff --git a/Runge_Kutta.m b/Runge_Kutta.m new file mode 100644 index 0000000..0106fe5 --- /dev/null +++ b/Runge_Kutta.m @@ -0,0 +1,23 @@ +function res=Runge_Kutta(fx,a,b,alpha,N) +h=(b-a)/N; +syms x y; +k1=h*fx; +k2=sliplify(h*subs(subs(fx,x,x+h/2),y,y+k1/2)); +k3=sliplify(h*subs(subs(fx,x,x+h/2),y,y+k2/2)); +k4=sliplify(h*subs(subs(fx,x,x+h),y,y+k3)); + +x_list=zeros(1,15); +y_list=zeros(1,15); +x_list(1)=a; +y_list(1)=b; + +coefficient=simplify(y_list(1)+1/6*(k1+2*k2+2*k3+k4)); +for i=2:N + x_list(i)=x_list(1)+(i-1)*h; + y_list(i)=subs(subs(coefficient,x,x_list(i)),y,y_list(i-1)); +end +disp([x_list;y_list]); +res= y_list; +end + + \ No newline at end of file From b8b4aee7d60d74031a6c0a7d784f4342f2bb0881 Mon Sep 17 00:00:00 2001 From: geniusjoe <843520313@qq.com> Date: Wed, 10 Apr 2019 20:41:28 +0800 Subject: [PATCH 14/19] add runge_kutta --- Romberg_calculus.m | 2 +- Runge_Kutta.m | 47 +++++++++++++++++++++++++++++++++++++--------- 2 files changed, 39 insertions(+), 10 deletions(-) diff --git a/Romberg_calculus.m b/Romberg_calculus.m index c5ed148..4155810 100644 --- a/Romberg_calculus.m +++ b/Romberg_calculus.m @@ -4,7 +4,7 @@ % \param b: right bound % \return res: the fx calculas % -% can change precision and max_iteration_times and min_diff_value at the following +% can change precision at the following eps=1e-6; buf=zeros(15); % size % ------------------------------------------------------------% diff --git a/Runge_Kutta.m b/Runge_Kutta.m index 0106fe5..b30de32 100644 --- a/Runge_Kutta.m +++ b/Runge_Kutta.m @@ -1,18 +1,47 @@ -function res=Runge_Kutta(fx,a,b,alpha,N) +function res=Runge_Kutta(fx,a,b,alpha) +% \param fx: syms expression +% \param a: left bound +% \param b: right bound +% \param alpha: y0 value +% \return res: y_value list +% +% can change iterate times at the following +N=10; +x_list=zeros(1,15); % size +y_list=zeros(1,15); % size +% ------------------------------------------------------------% +% example: +% syms x y; +% fx=y; +% a=0; +% b=1; +% alpha=1; +% Runge_Kutta(fx,a,b,alpha,N); +% +% return +% the answer is +% 0 0.1000 0.2000 0.3000 0.4000 0.5000 +% 1.0000 1.1052 1.2214 1.3499 1.4918 1.6487 + +% 0.6000 0.7000 0.8000 0.9000 1.0000 0 +% 1.8221 2.0138 2.2255 2.4596 2.7183 0 + +% ------------------------------------------------------------% + h=(b-a)/N; syms x y; + k1=h*fx; -k2=sliplify(h*subs(subs(fx,x,x+h/2),y,y+k1/2)); -k3=sliplify(h*subs(subs(fx,x,x+h/2),y,y+k2/2)); -k4=sliplify(h*subs(subs(fx,x,x+h),y,y+k3)); +k2=simplify(h*subs(subs(fx,x,x+h/2),y,y+k1/2)); +k3=simplify(h*subs(subs(fx,x,x+h/2),y,y+k2/2)); +k4=simplify(h*subs(subs(fx,x,x+h),y,y+k3)); -x_list=zeros(1,15); -y_list=zeros(1,15); x_list(1)=a; -y_list(1)=b; +y_list(1)=alpha; -coefficient=simplify(y_list(1)+1/6*(k1+2*k2+2*k3+k4)); -for i=2:N +coefficient=simplify(y+1/6*(k1+2*k2+2*k3+k4)); +disp([k1,k2,k3,k4,coefficient]); +for i=2:N+1 x_list(i)=x_list(1)+(i-1)*h; y_list(i)=subs(subs(coefficient,x,x_list(i)),y,y_list(i-1)); end From 09b36352e8d3cf48b88707bd6593cb386a2c90fd Mon Sep 17 00:00:00 2001 From: geniusjoe <843520313@qq.com> Date: Wed, 10 Apr 2019 20:46:20 +0800 Subject: [PATCH 15/19] Update README.md --- README.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/README.md b/README.md index 1cc7cbb..2780a39 100644 --- a/README.md +++ b/README.md @@ -10,12 +10,16 @@ |列主元高斯消去法| | `GE_equ.m` | |改进平方根法|对称矩阵的三角分解| `LDL_equ.m` | |Household方法|QR分解| `Household_QR.m`
`Household_equ.m` | +|龙贝格积分| | `Romberg_calculus.m` | + ## 迭代解法 |算法名称 | 文件名 | |:--------:|:--------:| |共轭梯度法 |`CG_equ.m`| |牛顿迭代法 |`newton_iteration.m`| +|四级RK法 |`Runge_Kutta.m`| + ## 插值 | 算法名称 | | 文件名 | From 5c2614ec5935b767932d82f6c59aef092c6f154b Mon Sep 17 00:00:00 2001 From: geniusjoe <843520313@qq.com> Date: Wed, 10 Apr 2019 20:57:25 +0800 Subject: [PATCH 16/19] Update Runge_Kutta.m --- Runge_Kutta.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Runge_Kutta.m b/Runge_Kutta.m index b30de32..36f19b4 100644 --- a/Runge_Kutta.m +++ b/Runge_Kutta.m @@ -16,7 +16,7 @@ % a=0; % b=1; % alpha=1; -% Runge_Kutta(fx,a,b,alpha,N); +% Runge_Kutta(fx,a,b,alpha); % % return % the answer is From de9f24c6d0812dfffa3cdd9dec143d58dd8f6808 Mon Sep 17 00:00:00 2001 From: geniusjoe <843520313@qq.com> Date: Thu, 11 Apr 2019 22:48:59 +0800 Subject: [PATCH 17/19] modify romberg (faster) --- Romberg_calculus.m | 46 ++++++++++++++++++---------------------------- Runge_Kutta.m | 2 +- 2 files changed, 19 insertions(+), 29 deletions(-) diff --git a/Romberg_calculus.m b/Romberg_calculus.m index 4155810..63803f9 100644 --- a/Romberg_calculus.m +++ b/Romberg_calculus.m @@ -2,8 +2,7 @@ % \param fx: syms expression % \param a: left bound % \param b: right bound -% \return res: the fx calculas -% +% \return res: the fx calculus % can change precision at the following eps=1e-6; buf=zeros(15); % size @@ -12,40 +11,31 @@ % fx=@(x)x^2*exp(x); % a=0; % b=0.8; -% Romberg_new(fx,a,b); +% Romberg_calculus(fx,a,b); % % return % the answer is % 0.5697 -% 2.0890 -% 0.2099 -% 0.3196 -% 0.3153 -% 0.3148 +% 0.3172 +% 0.3146 % 0.3146 % 0.3146 % ------------------------------------------------------------% i=1; -h=(b-a); -buf(1,1)=(b-a)./2.*(fx(b)-fx(a)); -buf(1,2)=h/2*(fx(b)+fx(a)+2*fx(a+h*i)); -buf(2,1)=(4*buf(1,2)-buf(1,1))/(4-1); -if (abs((buf(i+1,1)-buf(i,1)))=eps) - i=i+1; - n=2.^i; - h=(b-a)/(2^i); - buf(1,i+1)=h/2*(fx(a)+fx(b)); - for j=1:n-1 - buf(1,i+1)=buf(1,i+1)+h/2*(2*fx(a+j*h)); % get T_1_i+1 value - end - - for m=2:i+1 - k=i+2-m; - buf(m,k)=(4^m*buf(m-1,k+1)-buf(m-1,k))/(4^m-1); - end +buf(1,1)=(b-a)/2*(fx(b)+fx(a)); +while (1) + i=i+1; + buf(1,i)=1/2*buf(1,i-1); + for j=1:2^(i-2) + buf(1,i)=buf(1,i)+(b-a)/(2^(i-1))*fx(a+(2*j-1)*(b-a)/2^(i-1)); % get T_1_i+1 value + end + + for m=2:i + k=i+1-m; + buf(m,k)=(4^(m-1)*buf(m-1,k+1)-buf(m-1,k))/(4^(m-1)-1); + end + if (abs((buf(i,1)-buf(i-1,1))) Date: Thu, 11 Apr 2019 23:35:45 +0800 Subject: [PATCH 18/19] Update Romberg_calculus.m --- Romberg_calculus.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Romberg_calculus.m b/Romberg_calculus.m index 63803f9..74c5700 100644 --- a/Romberg_calculus.m +++ b/Romberg_calculus.m @@ -38,6 +38,6 @@ break; end end -disp(buf(:,1)); +fprintf('%.6f\n',buf(:,1)); res=buf(i,1); end \ No newline at end of file From ce6682cc6ef501391842c2fd6aba6a2d2aa35573 Mon Sep 17 00:00:00 2001 From: geniusjoe <843520313@qq.com> Date: Wed, 17 Apr 2019 20:41:41 +0800 Subject: [PATCH 19/19] Update Runge_Kutta.m --- Runge_Kutta.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Runge_Kutta.m b/Runge_Kutta.m index 70fded7..c9e4660 100644 --- a/Runge_Kutta.m +++ b/Runge_Kutta.m @@ -6,7 +6,7 @@ % \return res: y_value list % % can change iterate times at the following -N=5000; +N=20; x_list=zeros(1,15); % size y_list=zeros(1,15); % size % ------------------------------------------------------------% @@ -43,7 +43,7 @@ disp([k1,k2,k3,k4,coefficient]); for i=2:N+1 x_list(i)=x_list(1)+(i-1)*h; - y_list(i)=subs(subs(coefficient,x,x_list(i)),y,y_list(i-1)); + y_list(i)=subs(subs(coefficient,x,x_list(i-1)),y,y_list(i-1)); end disp([x_list;y_list]); res= y_list;