真實解 y(t)= sqrt(1+(1.0/exp(100*t)))
/* ex5-11.c based on Four-Order Runge-Kutta
* Method to approximate the solution of the
* initial-value problem
* y'=f(y,t), a<=t<=b, y(a)=y0
* at (n+1) equally spaced numbers in the interval
* [a,b]: input a,b,n,and initial condition y0.
*/
// y'(t)=(50.0/y)-50*y , 0<=t , y(0)= sqrt(2)
// 真實解 y(t)= sqrt(1+(1.0/exp(100*t)))
#include <stdio.h>
#include <math.h>
#define F(y,t) ((50.0/y)-50*y)
#define W(t) (sqrt(1+(1.0/exp(100*t))))
void main()
{
int i,n=10;
double a=0.0,b=1.0,y0=sqrt(2),k1,k2,k3,k4,h,t,y,err;
h=(b-a)/n;
t=a;
y=y0;
err=fabs(y-W(t));
printf("t y(t) w(t) error\n");
printf("=====================================\n");
printf("%.2lf %10.7lf %10.7lf %10.7lf\n" ,t,y,W(t),err);
for(i=1;i<=n;i++)
{
k1=h*F(y,t);
k2=h*F((y+k1/2.0),(t+h/2.0));
k3=h*F((y+k2/2.0),(t+h/2.0));
k4=h*F((y+k3),(t+h));
y=y+(k1+2*k2+2*k3+k4)/6.0;
t=a+i*h;
err=fabs(y-W(t));
if(i%1==0)
printf("%.2lf %10.7lf %10.7lf %10.7lf\n" ,t,y,W(t),err);
}
return;
}
輸出畫面
t y(t) w(t) error
==================================================
0.00 1.4142136 1.4142136 0.0000000
0.10 -15.8525966 1.0000227 16.8526193
0.20 -215.7458639 1.0000000 216.7458639
0.30 -2957.4012711 1.0000000 2958.4012711
0.40 -40541.0340389 1.000000 40542.0340389
0.50 -555750.0076712 1.0000000 555751.0076712
0.60 -7618406.3551152 1.0000000 7618407.3551152
0.70 -104435653.7847015 1.0000000 104435654.7847015
0.80 -1431638753.9652824 1.0000000 1431638754.9652824
0.90 -19625381252.2740898 1.0000000 19625381253.2740898
1.00 -269031267999.9240723 1.0000000 269031268000.9240723
Command exited with non-zero status 11
剛性方程（stiffness equation）是指一個微分方程，其數值分析的解只有在時間間隔很小時才會穩定，只要時間間隔略大，其解就會不穩定。
/* ex5-11.c based on Four-Order Runge-Kutta
* Method to approximate the solution of the
* initial-value problem
* y'=f(y,t), a<=t<=b, y(a)=y0
* at (n+1) equally spaced numbers in the interval
* [a,b]: input a,b,n,and initial condition y0.
*/
// y'(t)=(50.0/y)-50*y , 0<=t , y(0)= sqrt(2)
// 真實解 y(t)= sqrt(1+(1.0/exp(100*t)))
#include <stdio.h>
#include <math.h>
#define F(y,t) ((50.0/y)-50*y)
#define W(t) (sqrt(1+(1.0/exp(100*t))))
void main()
{
int i,n=1000;
double a=0.0,b=1.0,y0=sqrt(2),k1,k2,k3,k4,h,t,y,err;
h=(b-a)/n;
t=a;
y=y0;
err=fabs(y-W(t));
printf("t y(t) w(t) error\n");
printf("=====================================\n");
printf("%.2lf %10.7lf %10.7lf %10.7lf\n" ,t,y,W(t),err);
for(i=1;i<=n;i++)
{
k1=h*F(y,t);
k2=h*F((y+k1/2.0),(t+h/2.0));
k3=h*F((y+k2/2.0),(t+h/2.0));
k4=h*F((y+k3),(t+h));
y=y+(k1+2*k2+2*k3+k4)/6.0;
t=a+i*h;
err=fabs(y-W(t));
if(i%100==0)
printf("%.2lf %10.7lf %10.7lf %10.7lf\n" ,t,y,W(t),err);
}
return;
}
輸出畫面
t y(t) w(t) error
=====================================
0.00 1.4142136 1.4142136 0.0000000
0.10 1.0000227 1.0000227 0.0000000
0.20 1.0000000 1.0000000 0.0000000
0.30 1.0000000 1.0000000 0.0000000
0.40 1.0000000 1.0000000 0.0000000
0.50 1.0000000 1.0000000 0.0000000
0.60 1.0000000 1.0000000 0.0000000
0.70 1.0000000 1.0000000 0.0000000
0.80 1.0000000 1.0000000 0.0000000
0.90 1.0000000 1.0000000 0.0000000
1.00 1.0000000 1.0000000 0.0000000
Command exited with non-zero status 233
h=0.1
h=0.01
沒有留言:
張貼留言